{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from matplotlib.patches import Arc, Rectangle\n",
    "from matplotlib.lines import Line2D\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ls_circle(points):\n",
    "    '''\n",
    "    Input: Nx2 points\n",
    "    Output: cx, cy, r\n",
    "    '''\n",
    "    xs = points[:,0]\n",
    "    ys = points[:,1]\n",
    "\n",
    "    us = xs - np.mean(xs)\n",
    "    vs = ys - np.mean(ys)\n",
    "\n",
    "    Suu = np.sum(us**2)\n",
    "    Suv = np.sum(us*vs)\n",
    "    Svv = np.sum(vs**2)\n",
    "    Suuu = np.sum(us**3)\n",
    "    Suvv = np.sum(us*vs*vs)\n",
    "    Svvv = np.sum(vs**3)\n",
    "    Svuu = np.sum(vs*us*us)\n",
    "\n",
    "    A = np.array([\n",
    "        [Suu, Suv],\n",
    "        [Suv, Svv]\n",
    "    ])\n",
    "\n",
    "    b = np.array([1/2.*Suuu+1/2.*Suvv, 1/2.*Svvv+1/2.*Svuu])\n",
    "\n",
    "    cx, cy = np.linalg.solve(A, b)\n",
    "    r = np.sqrt(cx*cx+cy*cy+(Suu+Svv)/len(xs))\n",
    "\n",
    "    cx += np.mean(xs)\n",
    "    cy += np.mean(ys)\n",
    "\n",
    "    return np.array([cx, cy]), r\n",
    "\n",
    "\n",
    "def project_point_to_circle(point, c, r):\n",
    "    direction = point - c\n",
    "    closest = c + (direction / np.linalg.norm(direction)) * r\n",
    "\n",
    "    return closest\n",
    "\n",
    "\n",
    "def make_arc(c, r, start=0, end=2*np.pi):\n",
    "    theta = np.linspace(start, end, 100)\n",
    "    print(start, end)\n",
    "    x1 = r * np.cos(theta) + c[0]\n",
    "    x2 = r * np.sin(theta) + c[1]\n",
    "\n",
    "    return np.stack([x1, x2], 1)\n",
    "\n",
    "\n",
    "def get_angle_plot(line1, line2, offset=1, color='k', origin=[0,0]):\n",
    "    l1xy = line1.get_xydata()\n",
    "    l2xy = line2.get_xydata()\n",
    "\n",
    "    u = l1xy[1] - l1xy[0]\n",
    "    v = l2xy[1] - l2xy[0]\n",
    "\n",
    "    angle1 = abs(np.degrees(np.arctan2(u[1], u[0])))\n",
    "    angle2 = abs(np.degrees(np.arctan2(v[1], v[0])))\n",
    "\n",
    "    theta1 = min(angle1, angle2)\n",
    "    theta2 = max(angle1, angle2)\n",
    "\n",
    "    return Arc(\n",
    "            origin, offset, offset,\n",
    "            0, theta1, theta2, color=color, linewidth=1.0, zorder=1.0)\n",
    "\n",
    "\n",
    "def arrow(ax, p1, p2, text='', color='k', head_width=0.02):\n",
    "    line = Line2D([p1[0], p2[0]], [p1[1], p2[1]], color=color)\n",
    "\n",
    "    ax.add_line(line)\n",
    "    ax.arrow(\n",
    "            p2[0], p2[1],\n",
    "            (p2[0] - p1[0]) * 0.001, (p2[1] - p1[1]) * 0.001,\n",
    "            head_width=head_width, color=color)\n",
    "\n",
    "    if text:\n",
    "        ax.text(p2[0], p2[1], text)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = {'butter1': (252, 233,  79),\n",
    "        'butter2': (237, 212,   0),\n",
    "        'butter3': (196, 160,   0),\n",
    "        'orange1': (252, 175,  62),\n",
    "        'orange2': (245, 121,   0),\n",
    "        'orange3': (206,  92,   0),\n",
    "        'chocolate1': (233, 185, 110),\n",
    "        'chocolate2': (193, 125,  17),\n",
    "        'chocolate3': (143,  89,   2),\n",
    "        'chameleon1': (138, 226,  52),\n",
    "        'chameleon2': (115, 210,  22),\n",
    "        'chameleon3': ( 78, 154,   6),\n",
    "        'skyblue1': (114, 159, 207),\n",
    "        'skyblue2': ( 52, 101, 164),\n",
    "        'skyblue3': ( 32,  74, 135),\n",
    "        'plum1': (173, 127, 168),\n",
    "        'plum2': (117,  80, 123),\n",
    "        'plum3': ( 92,  53, 102),\n",
    "        'scarletred1': (239,  41,  41),\n",
    "        'scarletred2': (204,   0,   0),\n",
    "        'scarletred3': (164,   0,   0),\n",
    "        'aluminium1': (238, 238, 236),\n",
    "        'aluminium2': (211, 215, 207),\n",
    "        'aluminium3': (186, 189, 182),\n",
    "        'aluminium4': (136, 138, 133),\n",
    "        'aluminium5': ( 85,  87,  83),\n",
    "        'aluminium6': ( 46,  52,  54),\n",
    "        'indigo': (114,  33, 188),\n",
    "        'maroon': (103,   7,  72),\n",
    "        'turquoise': ( 64, 224, 208),\n",
    "        'green4': ( 0, 139, 0)}\n",
    "\n",
    "\n",
    "for k, (r, g, b) in colors.items():\n",
    "    colors[k] = (r / 255, g / 255, b / 255)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.283185307179586 8.230972752405258\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAALQCAYAAAC5V0ecAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAuIwAALiMBeKU/dgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xuc3GV99//XlU12s+GwCadkc4CN4ZRFZICkRsHTDWq1d2m12iradqwi7S1tlYH6612hLdb7Vm+Gemu0xaBOT9AW71tLW08cRG+RIEEGhIRTSCCHDackS8JuNsnm+v3xnYXZ2ZnNzu7M7On1fDzmsZnv9d3ruhIS9r3XXt/PFWKMSJIkSRqZGeM9AUmSJGkyMUBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVUwQEuSJElVMEBLkiRJVTBAS5IkSVWYOd4TmG5CCDOBFYXXucBioBnYD2wF7gPWAetijAfHa56SJEkqL8QYx3sO00IIYQFwCXApsGgEn7INuB5YE2PcUc+5SZIkaeQM0HUWQmgCLgeuAWaPoot9wNXAdTHG/lrOTZIkSdUzQNdRCGEpcBPw2hp0txa4OMa4qQZ9SZIkaZQM0HUSQugEbgPay7Wf1tzM2a2tnN7SwpwQ6ImRR/r6uL+3l0f376/UbRdwYYxxfZ2mLUmSpMMwQNdBYeX5LkrCcwAubmvjknnzOKu1teLnP9Dby5pdu7ixu5sy/3W6gPNciZYkSRofBugaK+x5vouSbRvLmptZ3d7OqjlzRtzX2p4eLuvqYuPQFem1wPnuiZYkSWo860DX3uWUhOeVra3c3tFRVXgGWDVnDrd3dLBy6Gr1KuATY5qlJEmSRsUV6BoqlKrbRFG1jWXNzdze0UFbU9Oo++3u7+eCzZtLV6L3AUstcSdJktRYHqRSW5dQFJ4DsLq9fUzhGaCtqYnV7e2886mnivdEzwY+AvzVmDqfIDK5fACaSP5ODrz2ZdOpfSP43GOB40j+yAdeALHwouT9wOtQ0au/6Fp/4XUQ2D+SOUiSpOnDAF0jhRMGLy2+dnFbW9XbNipZNWcO729r48bu7uLLl4YQPjuRTizM5PIzgJbCq7nwmlXyceDXM0tepTYAT4xg2GOBU8c69wp2kuxpH1Yml28BTgIOkJwq2Vf8MZtO+aMeSZKmCAN07ayg5ITBS+bNq+kAH503rzRALyY5Dvyemg5UhUwufzpwDK+E5lk17H6ke/TrGU5H+qDmHOC0So2ZXH4gTPeRbL8p+zJoS5I08Rmga2dF8ZvTmpuHLVU3Gme1tnJqczOPDd4LvYIaBOhMLt9EEgIHXtuy6VTFgtRF2khWgOthpHtfJkKAPtw3DgMr70cNd1Mml98H9AC9hY89L2x/cv/fX/Xejnjo0AqSb5gWF/raD2wF7gPWAesm0k8jJEmaqgzQtXNu8Zuzaxyei/stCdDnVrq3VCEkHwEcWfRxTuHXLSW37wGeH0G3fSMdfxQmwgr0SANprVbeZxdePL3h3rl3f/tv3r5j0/q3x0OHjhvmc9KFj9tCCNcDa3y4VJKk+jFA187i4jent5Tm0dpYPrTfxaUXMrn8LOBoktXOI4te1aT6kW7erucDdpNpBbq5VgMe6Ns349/+9x//2pZH1n0gxkPV9LsIuAbCnx197IJr33X5l7PHLVr2Yjadsl64JEk1ZICunUFBZ04Ile4bk9aSfmc0zZyTyeWX8EpgPoqiSiBjcMQI7xvJNo/RGukKdC/wAoMrbMAr1TgGfh2Kfj1jmFdT0b0NXYF+7N5b59/29//jin17d1fcT314sWXPzmf+7ObPfuQ9F/zOf782k5RW3AO8WHh1Az3ut5YkaXQM0LUzKEj21Km+dm9Jv61HzTsCSNVhqFqsQEeSP5eByhTFvz5QeB2s8BooJXdY2XSqi+SI85opVBOZychXt/eRVOxoZpQPUz74o/+75If/9LlP9x/Yf0y59tOamzm7tZXTW1qYEwI9MfJIXx/39/by6NDTKundu/u076751Of29ey56jVvevcWYH5Rc38mly8O1C8CrlZLkjQCBuja2Vr85pG++mwN3lDSb/PsI0ayT3k0RhqgXwQeJwnF+3il0kRfNp06UKe51V02nTpEFavr2XTqaeDpgfeFAD7w4GBL0cfZJa9WYMZj9946v1x4DiTlEC+ZN2/Yh1If6O1lza5d3NjdPSjx9x/Yf8wP/+lzn54956hPnrryrc8UNTUB8wqvATGTy+8Bdhe99hT+LCRJUoEnEdZICOEy4EsD709rbmbtsmU1H+e1GzcOeohw6Wve8Lfv+sQXv1PDISJJEO7OplP31rBfVfD2D/9ly625T/8kHuofVMllWXMzq9vbq6olvranh8u6ukpPraT1yLmPfuTa735yVsvsasPwIZIV6uez6dQjVX6uJElTkivQtbOu+M2j+/fzQG9vTUvZ5Xt7SytwsPCUs0Zy0EipCLwE7C18fIlCyTSg1xXHxvrB1//ijygpg7iytZWblyyp+hTLVXPmcHtHB+/dsoV7e3tfvt67d/dp//bFj//ae678229VOb0ZJKvUbu2QJKnAAF0764BtFB2msmbXLlbXMECv2bWr9NLWsy9833AB+gDJw2N7S14+QDZBhBAWANcUX1vW3Dyq8DygramJm5cs4YLNmwetRG/ZcO8Hnt5w7w9PXL5y9yi6HfKXr5xCqcToN2GSpKnMAF0jMcaDhRq8L4ehG7u7+eDcuTU5zvvunh5uGnwKIcD1zbOP2EOyl3YPr1Ra2EOyd7WeJeZUG5dQVDUlAKvb20cdnge0NTWxur2ddz711Mt7omM81Hz3t//mbScuX/mvo+hy5wjvWwKckcnld5FURnkB2OXDiZKkqcQAPUaFFbdjgON/7Y//9y9u+eIn9g/U7o3AZV1d3N7RMaZA1N3fzx92dZWWg9gH3ECyMrjfFeXJJ4QwE7i0+NrFbW01+YYLku0c729rG3T8+7bH7n/zjk3rP75gaecRJKdIHl14He7/BSNagSY5lXJG4ePACZWHMrn8TuC5wutF/75KkiYzHyIcpUwu3wYs55XAAMA3P3/pu57e8LMPFd872v2skITn0v2sBVfGGK+tfuaaKEIIq4C7i6/d2dFR033zD/T28ubNm0svr4oxvnz8eyaXDyR1v+cWvdp45e/1nmw6dedIxsvk8m9j6KmWpfaTnHL5HMnDiT0j6VuSpInCAD1KmVz+COC/lF4/0Ldvxpor3vG50oMwallRAVgLnB9j9Mfik9h4VW4BLosxfnm4zymU4TuKJEwfyqZTWw43TiaXPxJ4yyimuBd4FngG2On+aUnSRGeAHoNMLn8hZY7HfuzeW+d/d82nPleupu/729r46Ahq+n511y5uKqnpW9AFnBdj3DTW+Wt8hRC+AaQH3r+vrY2/Wbiw5uP8/vbt/Mvg/fPfiDH+Xq3HyeTyJwGvGWM3B0lWp58BnnUfvyRpInIPdEFhxe1Y4ARgwwhXwZ4FTiq9eOrKtz6zr2fPVaUHY0SSBwtv7O7m1MKpcstbWmgNgd4Y2VA4Va60VF2RLuBCw/OUsbj4zekth9v5MDrLh/a7uNx9NTDv8Lcc1kxgQeFF4bTELmBHNp16sQb9S5I0ZtM6QGdy+ZnA8UA7yTHHA38ez5Lszzyc5ygToAFe86Z3P36wr/d3/t83v3hN/4H955S2P7Z//3BBuZy1wMWG5ymlufjNnBDqMkjr0H7rk9ThAWATrzxAeCyjONK8xMBDjqdlcvkfZNOp+hzxKUlSFaZdgM7k8s0kYXkByWrzjDK3zWdkAfp5koXlQFJz+bnCteez6dRLpFOEm679HnA5SXm72RV7qmwfcBXw1+55nnIGfQfVU6ftVL1D+61LCC1U1uguvJ4sPJzYBhxH8o3qMZT/9zYSOw3PkqSJYloE6EJoXgAsJPlifrilvvnAQ4frN5tOHcjk8g+S1F7uLleaqxB6/1cI4R+Aj5CULRvJj9C3AtcDN8QYd4zgfk0+W4vfPNJXn3y4YWi/W8vdV2uFfw+7C68niks+Fl5HV9Gd/wYkSRPGlH2IMJPLz+KV0Hw8hw/NpX6YTaf21npehdq/55Ic3XwuSZhuIVkV3ArcR3Kq4X0xxoO1Hl8Tx0SuwtEImVx+Nsm/zfmFj8N9Q39HNp16aQR9zif5N78deM6KHpKkephSK9CFFa6B0Fxpe8ZIzScpr1VThVB8T+Gl6W1d8ZtH9+/ngd7emtaBzpd/KHVduXsbrVBhYwuwpfAQ7zySf3cnkJTQG7BnJOG5YDHJv//FwIFMLt9F8o3pTg9vkSTVyqQP0IV9lseRfMFsB8Z2BnJSRutZkm0ZUj2tA7YBiwYurNm1i9U1DNBrdg05QHDgpxwTSmGleODo7/WZXH4OrzyrMKJjxAvfQM8vujQLOLHw2pfJ5bcBW63mIUkaq0m7haNQQeNUkvAxmofzivWR7LHsAl7wx75qlBDCVSQPmCbvge+cdFJNjvO+u6eHX3nqqdJa4lfFGP9qzJ1PQJlcfiHJtqjD2UPyjcRW60xLkkZjMgfoAFzI6MNzD6+E5l3+eFfjIYSwgKT028t/j5c1N3N7R8eojn4f0N3fzwWbN5eeYrkPWDpVH0rN5PLnkmzfGKlIUjlnC0mdab9xliSNyKQN0ACZXL4TqOapqx6Sh4u2Z9Op7sPdLDVCCOFK4PPF11a2tnLzkiWjCtHd/f28d8sW7u3tLW26MsZ47ehnOnEVtm+8ndFv4TpAsp1mSzad2l2ziUmSpqTJHqCPBt50mNv2kYTmbX5h1EQUQmgC7gJeW3x9WXMzq9vbq9rOsbanh8u6ukpXniE5iOf8qVpLvPAQ4nEkK9ALGNsBLnuAp0i2eByowfQkSVPMpA7QAJlc/k0MrSe7n0Joxu0ZmgRCCEtJQnT7oOvA+9va+Oi8ecNW53igt5ev7trFTd3dpXueIdmmdN50OcWyEKbnkzxYPJZqPIdI/j/yVDadGtGDjJKk6WEqBOhlQCfJF7sdJA8HWf9Vk04IoRO4jZIQPeDU5mbObm1leUsLrSHQGyMb+vq4v3ypugFdwIUxxvV1mvaEVqgH307ysPFxY+hqQzadeqI2s5IkTXaTvowdSWDeT/IQkD9u1aQVY1wfQjgPuBFYVdr+2P79wwXlctYCF0+XledyCv9PeBp4unBwy2JgCXBklV111XpukqTJa9KvQEtTTWFP9OUk5e1GU2VmH3AV8NdTdc/zWGVy+XkkQXohh98v/Xw2nbq7/rOSJE0WBmhpgiqUuPsIcCnJyunhbAWuB26YqqXqaq3o9NIlJMeJl3NfNp3a3rhZSZImOgO0NMGFEGaSHBCyovBxMdBCcgDQwMmC64D7CkfFaxQyuXwrr5xcOLDyvx+4dSTPVBSqAp0MbLRMpiRNbQZoSSpSOKTpBOAk4MVsOvXICD8vRbKSDcmR5BuBZ60CJElTjwFaksao8IDiBQwtmfcS8CTJAS3uR5ekKcIALUljlMnlTwdOGeaWAyRHtm/KplNVlVKRJE08BmhJGoPCg4hvZWSnH/aTnHL4ZDadGnLWuiRpcjBAS9IYZHL5JUCqyk+LJA+AbsymU3tqPytJUj0ZoCVpDAoPHc4HlgHHjKKLHcDj2XRqd00nJkmqGwO0JNVIJpefC7yK5ICWUOWnPws8lk2ndg13U6Gs4QoGlzVsJim5V1zWcJ1lDSWpPgzQklRjhZrSS0lK4c2s8tOfIwnSO4svFg7WuYTkYJ1FI+hnG8nBOms8WEeSassALUl1ksnlZ5GE6FeRHH5TjeeAx6770NndjP1o96uB6zzaXZJqwwAtSXWWyeVnkByycjIwZ6Sft/6ufz/yezdc/cfAa2swjbXAxTHGTTXoS5KmNQO0JDVI4YHDhSRB+ujh7n3wR/93ye1/95mrYzw0v1z7ac3NnN3ayuktLcwJgZ4YeaSvj/t7e3l0f8VS013AhTHG9WP5fUjSdGeAlqRxkMnl5wOnAnNL2x6799b53/3qn322/+CBY4uvB+DitjYumTePs1pbK/b9QG8va3bt4sbubsr8H74LOM+VaEkaPQO0JI2jTC5/PHAaMA/gQN++GWuueMfn9u3dfVrxfcuam1nd3s6qOSPeAcLanh4u6+pi49AV6bXA+e6JlqTRMUBL0gRQCNKnfvPzl3746Q0/+1Bx28rWVm5esoS2pqaq++3u7+e9W7Zwb++Qgw+vjDFeO/oZqxqWH5SmFgO0JE0QhVJ1myiqtrGsuZnbOzpGFZ4HdPf3c8HmzaUr0fuApZa4qy/LD0pT04zxnoAk6WWXUBSeA7C6vX1M4RmgramJ1e3tpSe7zAY+MqaOVVEIoSmEcCXJN0TXMLLwTOG+a4BNIYQrQwhj+48vqS4M0JI0ARR+xH9p8bWL29qq2vM8nFVz5vD+trbSy5cWxlUNhRCWAncBn2d0tbspfN7ngZ8U+pM0gbiFQ5ImgBDCKuDu4mt3dnQMW22jWg/09vLmzZtLL6+KMd5Ts0GmuRBCJ3Ab0F6u3fKD0tTgyoMkTQwrit+c1txc0/AMcFZrK6c2N/PY4KC2AjBA10BhpXhIeK5B+cF24LYQguUHpQnCLRySNDGcW/zm7BqH52H6PbfcfapOYa/yTZSE52XNzXznpJNYvXDhYb8hOqu1ldULF/Kdk05iWXNzaXM7cKN7oqWJwQAtSRPD4uI3p7e01GWQ5UP7XVzuPlXtckqOXF/Z2srtHR1V72NfNWcOt3d0sHJo4F4FfGJMs5RUEwZoSZoYBi05zgmh0n1j0jq03/ok9WmkUKrumuJry5qbR127G5LKKTcvWVJuJfrThfEkjSMDtCRNDIM2JvfU6QHv3pJ+Z7W0zsjk8vVJ69OH5QelacYALUkTw9biN4/09dVlkA0l/R45b/4+4E2ZXP7Yugw4xVl+UJqeDNCSNDHcV/zm/qFHb9dEab9zT1jyBHAU8PpMLn92Jpd3S0d1VlBySMol8+bVdICPDu1vMT78KY0rA7QkTQzrit88un8/D9Q4ROd7e0tL2LHwlLOeKHq7GHhjJpf3a8PINaz84HDjSmos/ycpSRPDOmBb8YU1u3bVdIDS/mbOann+7Avf90TJbZuy6dShmg48tVl+UJqGDNCSNAHEGA8C1xdfu7G7m7U9PTXp/+6eHm7q7h50bf7SM77XPPuI4rC8B3iyJgNOH5YflKYhA7QkTRxrgH0DbyJwWVcX3f39Y+q0u7+fP+zqGny6XQh9r/v13/9Bya0PuvpcNcsPStOQAVqSJogY4w7g6uJrG/fv571btow6RHf39/PeLVvYWLL3mRg/deLyld8B9hauPJ1Np3aOapDpbVzKDwL1KdMiaUQM0JI0sVwH3FN84d7eXi7YvLnq7Rxre3q4YPNm7h36MOJa4K+z6dQLwI+ADcD6kfabyeUtofaKcSk/WDqupMYKsU7fLUuSRieEsBS4C2gfdB14f1sbH503b9hKDw/09vLVXbu4qbubMv+H7wLOizFuGs3cCoeuvIFk5fqhbDq1/zCfMqWFEC4DvjTw/rTmZtYuW1bzcV67cWNpBZXLYoxfrvlAkkbEAC1JE1AIoRO4jZIQPeDU5mbObm1leUsLrSHQGyMb+vq4v0ypuiJdwIUxxhGvNpfK5PInA8sLb/eT7JvuGm1/k10IYRVwd/G1Ozs6alrKLt/by1s2by69vCrGeE+Z2yU1gAFakiaowkr0jcCqGnS3Frh4tCvPAJlc/kjgTQzd/rcN+EU2nTowhvlNSoUTATdTdJjKB9raWL1wYc3G+Nj27dw4uILKVmBpoXKLpHHgHmhJmqAKYfd84E8oqs5RpX3AlcD5YwzPATiL8l83FgFvyeTy80fb/2Q1HuUHgesNz9L4cgVakiaBEMIC4CPApYysBvBWkmB3Q6G6x5hkcvkO4MwR3Po08HA2nZo2Aa/w32YTMHvg2rLmZm7v6KCtqWnU/Xb393PB5s2lFVT2kaw+j/m/qaTRM0BL0iRS2DJwLslRzueShOkWkrJmW4H7SE41vK+Wq5SZXH4ekAKOHMHtvUA+m049X6vxJ7p5TU1X7T506JriaytbW7l5yZJRheiB8oNlKqhcGWO8dgxTlVQDBmhJ0ohkcvkZwKnAySRFQQ5nI/DIVD+cZXdnZ2vfoUMfXvnkk//flgMHFhW3LWtuZnV7O6vmzBlxf2t7erisq2to7e5kH/v5McaxnawjacwM0JKkqmRy+bnA2YxsNfpF4OfZdGpPfWc1PnZ3ds4Cfhs48Z6enrnv2bLlw3sPHRr05zIRyg9Kqi0DtCSpaplcvgk4DRhJ0eNDJPuiN9d1Ug22u7NzBvCbwOkD127bu/f4D23b9julIXrAeJUflFRbBmhJ0qhlcvljSFajR7JH4VmSvdGT/hjq3Z2dAXgnsLK07Z6enrmXbN/+G1sOHBjJw56HM+byg5JqzwAtSRqTwtHencBJI7h9P3B/Np16tr6zqq/dnZ1vAC6o1N536NDzSx977KjeGK+mqDpHFfYBVwF/7Z7n6oQQFgK3VGi+KMa4vZHz0dRkgJYk1UQmlz+BpFJHywhufwJ4dDI+YLi7szMF/Powt+wBbpi7fn33eJcfnI4M0GoEA7QkqWYyuXwLyYErIzlUZRfJA4a1OXWkAXZ3dp4CvJ/KB5H1AV+fu379M8UXx6v84HRkgFYjzBzvCUiSpo7C/uafFQ5e6QSGK4I8D1hIsho94e3u7FwEvJfK4bkf+OfS8Awvn1h4T+ElaZLzKG9JUs0VKm78GBhyDnWR50lqRU94uzs7jwEuBpqHue1bc9ev92E/aRowQEuS6iKbTu0FfkL5kDzwMOGE30e4u7PzCOCDwBHD3Pa9uevXP9SgKUkaZ+6BliTVXeEBw7N5ZQX3nslQiWN3Z2czkCbZalLJT+euX/+DxsxIh+MeaDWCK9CSpLorhOUfAS8AT0yS8NxEsud5uPD8C+DWxsxI0kRhgJYkNUQ2ndoH3A08Ot5zOZzCQSm/CpwyzG2bgH+bu369P8qVphmrcEiSGqaw53nEgTOTy58I7M2mUzvrN6uy3kJS07qSHSQVNyw5J01DrkBLkiakwjHhrwFel8nlR3LKYU3s7uxcAbxxuFuAf5q7fv2kP5J8ijoIdFV4+Q2PasKHCCVJE07hQJY3MvgY7KeBX9Tz9MLdnZ2nA78FhAq39AJfm7t+/fP1moOkic8VaEnShJLJ5QNwDoPDM8CJwOszuXzp9ZrY3dl5IvAeKofng8CNhmdJrkBLkiaUTC5/OsM/vLcPWJdNp3aVaywcm72CwcdmN5PUni4+NnvdwLHZuzs7jwd+D2itMGYE/mXu+vWPVP0bkjTlGKAlSRNKJpc/i2S1eTiHSLZzPD1wIYSwALgEuBRYNIKhtgHXX3TUUTf+3eLFFwFtw9z7H3PXr183gj4lTQMGaEnShFN4aPDVHH6r4ZP/+tmPPLL10fsuB65h6LaPw2qCA7/Z1nbnF9rb724OodwXxR/PXb/+jmr7lTR1uQdakjThZNOpp0hqRg9b6eKxe2993bNPP/IA8HlGEZ4B+mHWTd3db12xcePv3dPTM7ek+X7gh6PpV9LU5Qq0JGnCKjwwuAKYV9r24I/+75If/tPnPt1/YP8x5T73mIVLWdBxBscuWsasllYO9PXywraN7Nj8MDu3byo73pEzZuz9xqJFf3/hkUc+BzxOUuu5v4a/JUlTgAFakjShZXL5GcCZFO2LfuzeW+d/d82nPjckPIfAGeddROqC32J+x/KKfT6zeQP52/+Fh++6BUq+Dh45Y8bev1u06DP/5cgjr5u7fv3+mv5mJE0JBmhJ0qSQyeWXAmcc6NvXtOaKd3xu397dpxW3z51/Im//8F+y6JThDhAcbNvjeb7/tT9n9zNPD7o+A352CF4fY3T1WdIQBmhJ0qSRyeWP/9fPfuR/bn30vg8XX29f9hredflqZs85quo+9/Xs4VvXXUbXxgdLm66MMV47hulqHIQQmoD5FZqf8Zsi1YIBWpI0aRRK1W2i6IHBufNP5OKr/3FU4XnAvp493HjNB0tXovcBS2OMO0bdsRouhLAQuKVC80Uxxu2NnI+mJqtwSJImk0sorrYRAm//8F+OKTwDzJ5zFG//8F9CGHQI4WzgI2PqWNKUZICWJE0KhRMGLy2+dsZ5F1W153k4i05JccZ5v1p6+dLCuJL0MgO0JGmyWEHJCYOpC36rpgOkLnhf6aXFJMeBS9LLDNCSpMliRfGbYxYuHbZU3WjM71jOMe1Lhx1XkgzQkqTJYtBK8IKOM+oyyPylncOOK0kGaEnSZLG4+M2xi5bVZZDjFp087LiSZICWJE0WzcVvZrW01mWQmc2zSy+11GUgSZOWAVqSNFkMOlb7QF9vXQY5uH9f6aW+ugwkadIyQEuSJoutxW9e2LaxLoM8v+2JYceVJAO0JGmyuK/4zY7ND9dlkGc2rR92XEkyQEuSJot1xW92bt/EM5s31HSAZzavZ2fXpmHHlSQDtCRpslgHbCu+kL/9X2o6QJn+tuIKtKQSBmhJ0qQQYzwIXF987eG7bmHb4/ma9L/tsft5+K5/L718fWFcSXqZAVqSNJmsAV4pkxEj3//an7OvZ8+YOt3Xs4fvf/0vIMZBl4EbxtSxpCnJAC1JmjRijDuAq4uv7X7mab513WWjDtH7evbwresuY/czTw+6vui0c/7+8m/cv3fUk5U0ZRmgJUmTzXXAPcUXujY+yI3XfLDq7RzbHs9z4zUfpGvjg4Outx4599F3f+Im9Ve4AAAgAElEQVTL3wFen8nljxrrhCVNLSEO/nGVJEkTXghhKXAX0F7SwBnn/SqpC97H/I7lFT//mc0byN/+z8me55Kvg02zmne+45K/+uSpK9/6TOHSfuDubDr1Yk1/E6qLEMJC4JYKzRfFGLc3cj6amgzQkqRJKYTQCdxGaYguOKZ9KfOXdnLcopOZ2Tybg/v38fy2J3hmU9lSdUASnt/ygU9e9Zo3vXtLSdMBkhDdXdPfhGrOAK1GMEBLkiatwkr0jcCqsfbV0nrkw2/90NVfKFp5LnUAWJtNp3aPdSzVjwFajeAeaEnSpBVj3AScD/wJxdU5qrMPuPLst1189qkr3zrc8YazgFWZXL5tlONImiIM0JKkSS3G2B9j/F/AUuAqksNPRmJr4f6lMcZr7/729QeAnwHPDfM5s4DXZXL5o8cyZ0mTm1s4JElTSghhJnAusKLwcTHQAvTxysmC64D7yh2SksnlZwArgROGGWY/8NNsOjW2AtSqObdwqBEM0JIkFQkhNF/+jfsPkoTvBcPc2kcSoq0VLU0zbuGQJE17IYQZIYR3hhDuAG7IplOHSFaqKz1QCMmq9uszufyRDZnkNBNCeHMIIYYQ7hzvuUilDNCSpGkrhDA7hPBh4CHgP4G3AO8PISwuhOh1wLPDdNFCsid6Tv1nK2miMEBLkqalEMKpwGbgBqD41JWZwB8BFEL0vQz/YOFsYFF9ZilpIjJAS5Kmq41Apf3Ll4YQjoYRhejHsunU43WY37RQ2DrzvRDCySO8f0YI4Q8KnxPqPT+pHAO0JGlaijH2A9dVaD4a+PDAm2w61U8Sol8ouW99Np16tD4znDb+BHg78FAI4ZoQQmulG0MIv0RSavArhc95S2OmKA1mgJYkTWc5hobiAR8PIcwaeFMI0fcAOwuXHsymUxvrO71p4Z3Ap0iqmlwFPBxC+NXiG0IIx4YQvgqsJamO8mPgdTHGOxo9WQksYydJmuZCCNeQBLdyLo4x3lR8IZPLzwKOyaZTw1XoUJVCCMeRBOk/AJqBh4EzgB0kB9gcS/Kw55/GGP9jvOYpgQFakjTNhRDmA0+RVNQo9XNgRfSLZcOEEJYCfwW8HxjY47wF+HPg72KMh8ZrbtIAt3BIkqa1GOMzwD9UaD4HeHPjZiMgFl6l1/rHYS5SWQZoSZIqP0wIcMVYOs7k8hUfitMrQgjHhBCywCPAB4D1haYdwBHA3wE/DyG8/TD9LAwhrKvwWljX34SmDQO0JGnaizFuACrtq31nCOGM0fSbyeXbgf+SyeVPHPXkprjCYTafJCkreDnQBVwEXFa45VHgNJJ63a8BvhdCuC2EcM54zFcCA7QkSQOuHabt8mo7y+TyJwErSL7WviaTyy8Y7cSmuO8CnwVaSfY+d8YY/734hhjjCzHGS4DXkRyxfgGwLoRgGTuNCwO0JEmJH5Mc3V3OB0MIIw7AmVz+NJLV0gEBODeTyx8zhvlNVdcC3wfOjDFeFWPsrXRjjPEe4JeAjwG3AT9qzBSlwQzQkiQBhUoblVahm3llS8GwMrn8mcCpZZpmAL+UyeWPGN0Mp6YY43/GGH85xjii0xxjjIdijF+JMb7NihwaLwZoSZJe8X9IStqV899CCCMJv5WOB4eknvGqTC5frmSepEnCAC1JUkGM8SDwhQrN84APHa6PbDq1CRhuNXUOsDKTyzdVP0NJE4EHqUiSVCSEcBTJwR1tZZo3AafEGA9bkziTy58FDFd9YwewLptOTegvxCGE44GzgUUkfyZHAvtIVtq3AU8Cjxa++Rh3hVJ1t1RovijGuL2R89HUNHO8JyBJ0kQSY9wTQvhb4JNlmpcC7wK+OYKufkFSWeL4Cu0LgE6SI6snlBDCKSSr7e8j+T0fTm8I4T6SB/u+HWN8oJ7zk8abK9CSJJUIISwiWW2eVab5Z8CqkRzvncnlZwLnAUcPc9tDhW0f4y6E0EZSSu4PgNFuMdkWY1xcu1lVxxVoNYJ7oCVJKhFj3AbcWKH5l0hC8WFl06mDwD0kWx4qOSOTy59Q3QxrL4RwInAvSbWRsezPrlQKUJoyXIGWJKmMEMKZwIMVmr8dY3zXSPvK5PJHk4TuSlsnDwI/yaZTe6qbZW2EEOYCeeCkkqZDJAeX/ITkwcjdwAHgKJKtKWcA5xY+DvhUjPEz9Z5zJa5AqxEM0JIkVRBC+D7wtqJLLwFfA74QY6xq20Umlz8eeC3JoSrl9AD/L5tO7R/NXMcihPB1hlYYuQm4Ksa4cQSf3wH8Jsnq9SUxxu/Xeo4jZYBWIxigJUmqIITwVuAHJBUzvghcH2PcOdr+Mrn8EiA1zC07gbuz6VTDDggJIbQDWxm8rfMzMcZPjaKvmcCMGGPDvwkomoMBWnXnHmhJkiq7jWRltSPG+D/HEp4BsunUFmC4Fd1jGHwEeCP8BoPzwA7gL0fTUYzx4HiGZ6lRDNCSJFUQEzfHGPtq2O0GkpBayZJMLr+shuMdzvKS9w/GGA80cHxp0jFAS5LUQIWDU+4HXhzmtuWFPdONML/k/SkhBE9JlIZhgJYkqcEK5e1+BlRa2Q7AuZlc/ogGTGdvyfulwN+EEMrVwJaEAVqSpHGRTad6SeouV3pgcBawsAFTuaPMtUuAB0MIvxtCaESIlyYVA7QkSeMkm07tIqm/XCoCD2bTqccbMI2bSI4dL3U6kAOeDyF8J4TwsUK5OmnaM0BLkjSOsunUNgZX5ugDfppNp55qxPiFBwbfCfy8wi2zgXcAq4FNIYSHQgj/I4TQ2Yj5SRORAVqSpFEKiQtDCN8LIbx2DF1tAJ4jOenvx9l0akzl8qoVY9wKrAL+ENh8mNvPAP4UeDiEcGcI4Q11np404RigJUmqUgihOYTw2yTbL24F3g5kRttfoTLHfSQrz/tqM8vqxBgPxBhXA68C3gxcCzxAsp2kkjcBPwohVH3oijSZeRKhJElVCCG8F/hrYFFJ0yHglBjjk42fVf2EEI4BLgB+Dfh1oNJDhe+IMX6vYROroFCCr7Q034BnYoz9jZyPpiYDtCRJVQgh/ArwHxWavxRj/KNGzqeRQghzgcuBTwLNJc23xBh/rfGzkhrPAC1JUhVCCDOAhxh6gh9AD7BkrEd+T3QhhLcC32PwVtAnY4yNPEFRGjfugZYkqQoxxkNAtkLzHOD36z2HTC6/MJPLt9V7nEpijLcCPym5fHA85iKNB1egJUmqUghhNkm1inJ7bZ8BTooxVjplcNQyufwMkpXvV5Gsdv84m04dqPU4IxFCuJukcseAO2KMF4zHXKRGcwVakqQqxRj3AV+q0Dwf+ECtx8zk8rOB15GEZ0hWu8+q9TgjEUI4CVhZcrl0RVqaslyBliRpFEIIxwJPkwTZUuuBV8cafZHN5PLHAecALWWaH6z20JUQwseBrcC3qq1KEUI4imT/8+uLLkfg9BjjY9X0JU1WrkBLkjQKMcYXgK9XaO4EfrmGw51I+fAM8OpMLn90lf2lgZuBjSGEz4QQzhnJJ4UQfhm4l8HhGeBvDM+aTlyBliRplEIIy4DHgVCmuWZ7gjO5/EzgjVSuwfwSyX7owz7IF0JoBvYCs0qauoB1JEd67wB2kiy0HUty+uDbeWX7SLH7gTfGGPce/nciTQ0GaEmSxiCE8E3gNyo0nxtj/HktximsMr+Byj893pZNpw47VmG1+b5azAm4C/iVGGN3jfqTJgW3cEiSNDbXDtM26uO9S2XTqReBh4e5ZVEml18ygq4WA2Ot3PEScAXw5okWnkMIJ4QQ/r3C64Txnp+mBlegJUkaoxDCT4DzyjT1A6+KMT5dq7Eyufy5wMIKzf3Aj7Lp1EvD9VF4EPCCwmsVSTWP0i0dpQ4BDwB/D/xDYQ/4hBNCWAjcUqH5ohjj9kbOR1PTzPGegCRJU8C1lA/QTcAfU8OVaOBBYC7lq380Aedkcvm7sunUoUodxBj3AN8uvAghtAAnA8tIwvlRJA8t7gF2k9S8vs99zlLCFWhJksYohNAEbABOKdO8h+R475ptdcjk8nNJAnulrZiPZ9OpR2o13mTiCrQawT3QkiSNUaGW8nUVmo8CLqnleNl0ajdJYK/klEwuf2wtx5T0CgO0JEm18ffA8xXaPl4oH1dLm4Dnhmk/J5PLH25fs6RRMEBLklQDMcYe4CsVmhcBv1nL8bLpVATywP4Kt8wGXlPLMSUlDNCSJNXOl4G+Cm1XhBDKHbgyatl0ah9JZYxKFmZy+UW1HFOSAVqSpJqJMT4L/F2F5rNIysbVVDad2gE8NcwtZ2Zy+dm1HleazgzQkiTVVrmHCfcB1wNP1mnMh0mO5y5nFpCq07jStGSAliSphmKMj/JKGbXngb8ATowx/n6MsS4BOptO9QM/ByrVpj0+k8t31GNsaToyQEuSVHufAS4lCc5/GWMcrlpGTWTTqW7g0WFu6czk8kfUex7SdOBJhJIk1ViM8WfAz8Zh6CeABSQnFZY6QHK64LDHfEs6PFegJUmaIgql7e4HSo/x3gbcmU2ndjZ+VtLU41HekiRNMZlcfinwapIa0Q9m06mucZ5Sw3iUtxrBLRySJE09m4FmYHM2napUl1rSKLkCLUmSpgxXoNUI7oGWJEmSqmCAliSpwULijSGEbK2P95ZUfwZoSZIaJIQwM4TwW8A9wI+Ay4E3jO+sJFXLAC1JUp0VVpz/EHgc+GdgZVHzFeMzq0Qml5+ZyeVnjeccpMnGAC1JUp3F5In9dwAdZZp/NYRwemNnlMjk8vOBNwNnjsf40mRlgJYkqTGyw7R9omGzADK5/OxMLn8u8EtAK7Aok8sf38g5SJOZZewkSWqAwsOCPwdSZZr7gBNjjM/Wcw6ZXD4AJwHLGXoWRA/JaYX99ZxDvYUQmqm8ov6LGOP+Rs5HU5MHqUiS1AAxxhhCuBb4xzLNLcDHgD+v8zQCyTaScl//5wCnAhvqPIe6KgTk+8Z7HpraXIGWJKlBQgizgCeBxWWaXyBZhe6p5xwyufwxwHkVmiPw42w69WI95yBNdu6BliSpQWKMB4AvVGg+Fvjdes8hm07tBJ6q0ByA1xS2ekiqwAAtSVJjrQEqrfBeHkJoasAcNpDsuy5nHuVXyCUVGKAlSWqgGOOLwFcrNJ8MXFTvOWTTqQPAQ8Pc0mltaKkyA7QkSY33ReBghbaGHKySTae2A5WqfjQD41KbWpoMDNCSJDVYjHEL8C8Vml8fQnh9g6byEHCoQltHJpdva9A8pEnFAC1J0vgY7mCVTEMmkE69BDwxzC1n+kChNJQBWpKkcRBjvB+4vULzu0IIJzdoKk8AvRXa5gFLGjQPadKwDrQkSeMkhPDLwHcrNH8lxvixRswjk8svAFZWaN4P3FF48HDCCyEcC/zPCs1/GmN8oZHz0dTkCrQkSePn+1SuhvGhQhisu2w6tYPhHyg8pRHzqJEW4JwKr5ZxnJemEAO0JEnjJCY/Bq60F7oV+IMGTme4BwqXZnL5Ixo4F2lCM0BLkjS+bgK6KrT9YQhhdiMmUXigcGOF5hnAGY2YhzQZGKAlSRpHMcY+4EsVmk8APtjA6TwB7KvQNj+Tyx/fwLlIE5YBWpKk8Xc98FKFtkwIoSFfr7Pp1EHgkWFuOcOydpIBWpKkcRdj3Al8rUzTfuCnQCP3H28FdldoOwo4qYFzkSYkA7QkSRPDF3jlIb5dwP8AOmKMH44x7mnUJLLpVKRyZRCAV7kKrenOAC1J0gQQY9xEEqL/CDgxxvhnMcZKDxfWVTad2gVsK7l8CHgS+EkhZEvT1szxnoAkSUrEGBtyhPcIrQcWAE3AdmBDNp3qGd8pSRODJxFKkqSyMrn8EmBPNp2qtCd6wgkhLARuqdB8UYxxeyPno6nJAC1JkqYMA7QawT3QkiRJUhUM0JIkSVIVDNCSJElSFazCIUnSJBJCeB3QEmO8c7znMh5CCIHkYJmjSQ52Obrk18cVXgAB6CY5kEaqGQO0JEkTXAihCbgIuAJ4PfCLEMJZcYJUAsjk8scCO8dSHzqE0AosKnotLvl1OzCXJChXc5DLExigVWMGaEmSJqgQwhzgd4HLgZOLms4E3gr8YDzmNSCTy88DlgPHAj9n6OErg4QQZpAcBb685HVaoY968NRE1Zxl7CRJmqBCCKcCj1A+BN4aY3xbg6cEQCaXPxo4HZhfdPkl4M5sOnUIIIRwPLACOAfoJAnKpwOtjZ0tG4GBOtaWsVNNGKAlSZrAQgjfAn69QnMqxvhAI+eTyeWXAq8uvrZ393OzfnHn/zn5oZ/c0rbnha5TgZUkK82N0Ae8COwpvA4BpwIDAWcbsLfwawO0asItHJIkTWzXUjlAZ4DfaeBcAJ7rfm5r833f/6cztz+RP2vPCzuW977UvYwYa5Up9pCE3m3A1qJfbwNe4JWw/CKwJ8Y4aH/zYQ5SkWrCFWhJkiawQtWJnwKryjQfBJbGGLfWeQ4zSLZivBV4G3A+Y1uE2w1sANYXPm4ANgHbYowvjnGunkSounMFWpKkCSzGGEMI1wLfLNM8E/hD4JO1HjeEsIQkLL8NuIDRPeTXQ/Jw4QMMDsvPTJQKItJouAItSdIEVyhj9xjwqjLNLwJLxrpyWxhnKfAe4L0k+5irEVvmHPXUkfNOeIQY73hh+5PfAzbEGA+OdV7VKFQuuaBC8+0xxp5GzkdTkwFakqRJIITwMWB1hebLY4x/Pcp+TyYJze8Bzq3iUw+2HjXv0XkLTnqgfdmZG179hl9/7NiFr+ottPUCdwxU5JCmGgO0JEmTQAjhCOBp4JgyzU8DJ8cYD4ywr2XAb5GsNKeqmMajJLWnbwXuvPwb93cAHRXufTCbTj1VRd/SpDFjvCcgSZIOL8b4EvCVCs0nkqwgVxRCaA0hfCCE8EOS0/k+w+HD8y7gX4FLgI4Y4+kxxj+KMf57jHEP8DhJ2bhyTsnk8uYMTUn+xZYkafJYTeVjqa8oVOwYJIRwVgjhS8B24B+BNx9mjOeBNcDbgfkxxt+KMd4QYxyympxNp/aRrH6X00pyBLc05RigJUmaJGKMzwD/UKH5HOBNACGEo0MIl4YQ7gXywGXA3GG6fhb4W+BCoD3G+NEY4w9GuCXkCSqvQp+cyeU9SltTjnugJUmaREIIncDDFZrvJDn6+3eAOYfpaifwz8DNwP+LMfaPdk6ZXP5MKu+F/nk2ndo22r6licg60JIkTSIxxvUhhP8EfqVM85s5/BaNW4GvAd+OMfbVaFpPkBzdXW61+WSSUwSlKcMVaEmSJpkQwpuBH1bxKduArwPfiDFuqsecMrn82VTe83xPNp16th7jSuPBPdCSJE0+PwLuO8w9/cC3SFaqT4oxXl2v8FzwxDBtp9RxXKnhXIGWJGkSCiG8D7ipQvNPgd+MMTZ060Qml18JLKjQfFc2ndpZ7zmEEOYCH6/Q/IUY4+56z0FTn3ugJUmanL4JfI6kBnSp5cB4BMXHqRygTwHuacAc5gD/tULbVxmfPxdNMW7hkCRpEooxHgQqHd89D0g3bjaJbDq1m6SOdDknZHL5Ixs5H6leXIGWJGny+hrwF0BbyfWDjN8hJo8Dx1VoWwY8EEKYCawovM4lmWszySExW0n2d68D1hW+UZAmFPdAS5I0iYUQPgt8svB2D3A98MUY45bxmlMml38jQ0M9j6+7vf/fv3zFucClwKIRdLWN5PezJsa4YyRjhxAWArdUaL4oxrh9JP1Iw3EFWpKkye1LwHuArwA3xBhfHOf5AGwkORkRIPb17N1xw5W/8o6+nhf/OzC7in4WAdcA/z2EcDVw3VgOfJFqxQAtSdIkFmPcFkI4JU6sHyl3kayGP/cPf/6+/ueefjQHvHYM/c0GPg+8O4RwcZ3L8UmH5RYOSZJUc5lcPlz3obOXA7cB7eXuOa25mbNbWzm9pYU5IdATI4/09XF/by+P7t9fqesu4MIY4/pyjW7hUCMYoCVJUs2FEJYCd1ESngNwcVsbl8ybx1mtrRU//4HeXtbs2sWN3d2USSpdwHnlVqIN0GoEA7QkSaqpEEITSXgetG1jWXMzq9vbWTVnzoj7WtvTw2VdXWwcuiK9Fji/dE+0AVqNYB1oSZJUa5dTEp5XtrZye0dHVeEZYNWcOdze0cHKoavVq4BPjGmW0igZoCVJUs2EEBaQVM542bLmZm5esoS2pqZR9dnW1MTNS5awrLm5tOnThfGkhjJAS5I0DYQQzgkhVN50XDuXUFSqLgCr29tHHZ4HtDU1sbq9nTD48mzgI2PqWBoFA7QkSVNUCGFGCOGdIYQ7SE73++06jzeT5JCUl13c1lb1to1KVs2Zw/vbhpzPcmlhXKlhDNCSJE0xIYTZIYQPAw8B/wm8pdCUCSHU82v/CkpOGLxk3ryaDvDRof0tJjkOXGoYA7QkSVPP1cANwPKS66cC/7WO464ofnNac/OwpepG46zWVk4duhd6Rbl7pXoxQEuSNPV8Fah05PUVdRx30Erw2TUOz8P06wq0GsoALUnSFBNj3AzcXKH5DSGEsRyrPZzFxW9Ob2mpyyDLh/a7uNx9Ur0YoCVJmpqyw7Rl6jTmoL0Vc0KodN+YtA7ttz5JXarAAC1J0hQUY1wH3Fmh+TdCCK+qw7CDjgvsqdNpx71D++2ry0BSBQZoSZKmrmsrXJ8BfLwO420tfvNIX31y7Yah/W4td59ULwZoSZKmru8CGyq0/V4I4Zgaj3df8Zv7e3tr3H3Ffu8rd59ULwZoSZKmqBjjIeC6Cs1HUHLoSQ2sK37z6P79PFDjEJ3v7eWx/ftLL68rd69ULyHWaX+SJEkafyGE2cBTwAllmncAHTHGmuy1KJwIuJmiw1Q+0NbG6oULa9E9AB/bvp0bu7uLL20FlsYYDxbmcBTw/gqfflOMcU/NJqNpywAtSdIUF0L4FPDpCs2/F2P8Rg3Hugq45uX3wHdOOqkmx3nf3dPDrzz1FCXJ5aoY41+NuXOpCgZoSZKmuBDCscAWoNzJJuuBV8caBYIQwgJgEzB74Nqy5mZu7+igralp1P129/dzwebNbBy8fWMfyerzjlF3LI2Ce6AlSZriYowvAF+v0NwJ/HINx9pBcpT4yzbu3897t2yhu7/S4YjD6+7v571btpSGZ0hWnw3PajhXoCVJmgZCCMuAx0l2VZS6I8Z4QQ3HagLuAgadeLisuZnV7e1VbedY29PDZV1d5cLzWuD8GOPoUrk0BgZoSZKmiRDCN4HfqNB8Tozx/hqOtZQkRLcPug68v62Nj86bx1mt5XaUJB7o7eWru3ZxU3d36Z5ngC7gvBjjplrNV6qGAVqSpGkihLAKuLtC8z/FGD9Y4/E6gdsoCdEDTm1u5uzWVpa3tNAaAr0xsqGvj/vLl6ob0AVcGGNcX8u5StUwQEuSNI2EEO4CXl+mqZ/kgbwtNR5vKXAjsKoG3a0FLnblWePNhwglSZpeKh3v3QT8ca0HK4Td84E/IamaMRr7gCtJ9jwbnjXuXIGWJGkaKTzg9whwcpnmPcCSGGN3mbZajL0A+AjJCYiLD3v/jKaueKj/K8ANVtvQRGKAliRpmgkh/AHwlTJNh4B3xxj/rc7jzwTOBVYA5845+pjXAC1hRtPB5tlHPD/3hCVPLDzlrCfOest71335Y2/8UZV9exKh6s4ALUnSNBNCmAM8DRxbuPQS8DXgC+OxRSKTy58GnFqh+Y5sOvXSSPsKISwEbqnQfFGMcXu185NKuQdakqRpJsbYA3wZ2AH8Kcm2jT8ex/3F2ypcfxFoaeREpJFwBVqSpGkohHAkcCDG2DfecwHI5PJvAOaSPDC4FdiWTaderLYfV6DVCDPHewKSJKnxYox7x3sOJTYAEdiZTadc3dOE5gq0JEmaMlyBViO4B1qSJEmqggFakiRJqoIBWpIkSaqCAVqSJFUUQjhuvOcgTTQGaEmSNEhIvDWE8H3gsULJO0kFBmhJkgRACKE5hPDbQB74AfA2YB7we+M6MWmCsQ60JEkihBCAdcCZZZo/EUL4SozxYIOnRSaXD8AxQDtwdDad+mmj5yCVMkBLkiRijDGE8B+UD9AdwLuBf23EXDK5/AzgBGABMB9oLmo7ejQnFEq15BYOSZI04EvAgQptVxRWqRuhCVgBLKEoPBe0N2gOUkUGaEmSBECMsQv4xwrNK4HzGzGPbDp1AHihQvOCRsxBGo4BWpIkFbtumLYrGjYL2FHh+tGZXH5OA+chDWGAliRJL4sxPgR8r0LzRSGE0xo0lUoBGlyF1jgzQEuSpFLXDtN2eSMmkE2neoHuCs3ug9a4MkBLkqRSd5DUgi7nd0MIJzRoHl0Vrs/L5PKlDxdKDWOAliRJg8QYI5VXoVuAjzVoKpW2cQTg+AbNQRrCAC1Jksr5V2Brhbb/FkKo+4N82XRqD9BToXl+vceXKjFAS5KkIWKMB4AvVGg+DvidBk3lmQrXTyicUig1nCcRSpKkStYAVwNHl2m7PISwJsbYX+c5PAssLXN9FjAP2FlyvQf4jwp9VVrNlqpigJYkSWXFGF8MIawBMmWaTwF+Ffh2nafxPNBPcjphqfmUBOgY427gL+o8J01zbuGQJEnD+d/AwQptdT9YJZtOHSIJ0eU0qhqINIgBWpIkVRRj3AL8S4Xm80IIr2vANCrtgz46k8vPbsD40iAGaEmSdDjZYdrKbe+otWeHabOcnRrOAC1JkoYVY7wfuL1C87tCCMvqOX7hVMI9FZoN0Go4A7QkSRqJSgerzAA+3oDxn6tw/XjL2anRDNCSJGkkvg88VKEtHUI4os7jVwrQzZQvsyfVjQFakiQdVuF479K90PuAvwXOjTG+VOcpvAAcqtDmNg41lAFakiSN1E3ADpLV4D8HTowx/kGM8bF6D5xNp/oZemjKAAO0GsqDVCRJ0ojEGPtCCO8AHo0x9o7DFJ4jOUa81DGZXL4pm071hxDmABdU+PzbY4yeRqgxM0BLkqQRizHmx3H454DlJdd6SA5amUlyYvxLvYoAABAoSURBVOFcktXxcu7D47xVAwZoSZI0WbwI7CYpafc88EKhxJ3UUCF5JkCSJGnyCyEsBG6p0HxRjHF7I+ejqcmHCCVJkqQqGKAlSZKkKhigJUnS/9/evYZYWt8HHP/9z5yZvc6uujuia22MLYHaSTUgAQNt0VAbqBVqA5ZsU7zguoFapbNB6av2RfsmcaU1tCoxXdItRHpJ2WL7xsJCL6GlNgv1EmKlAcHLXrqrXd3V3Tn/vpip2eh51vntPOecOTOfDwh6fjPP8xf2xdfH/3n+rSulaAxWLX+4AYBWlAU/X0r524j4yqjXA4MioAGAZSmldEspt0fEv0XEwYi4JSJ2lVIuGunCYEAENABwwUopsxHxXxHxrYi4/pzR5ojYNZJFwYAJaABgOV6OiE0Ns/tLKVPDXAwMg4NUAIALVms9VUr5WkT8bp/xjoj4tYj45iDXMLfv0FREXBIRl3z2N35n5h+++QeDvB0IaABg2f44Ih6KiPV9ZntKKX9WWz65bW7foS0RcXUshPP7T8C3bNuxpc37QD+2cAAAy1JrPRIR+xrGn4yIXxjAbbsRcWV8YPtI6XS6pdPRNwyUP2AAQBseiYimp8xzA7jfmxHR6zfodLr+DzsDJaABgGWrtX4/Ig40jG8upfxMm/d7+I7r5mMhoj+kMzEx0ea94IMENADQlq+eZzaIp9D/0+/DTmfCE2gGSkADAG3554j414bZF0opP9by/Y73+3BhD3QpLd8L3iegAYBWLL5po+kpdDci7mv5lieaBrZxMEgCGgBo07cj4r8bZrtLKa29Zu7hO647FRHv9pt1OgKawRHQAEBraq3zEbG3YbwlIu5u+ZZ9n0IXAc0ACWgAoG1/Gg37kyPigVJKm1/y63ufzoQvEjI4AhoAaFWt9e1YOJ2wnx+PiM+3eLv+T6AXtXgfeJ//OgMABuFrEfHliJjqM/tyKeWplo73/pGAXr9x+szGrdtfjojo9XpH3jt18tQ54777pSGrtHw0PQBARESUUp6MiLsaxjfWWg+2cZ+5fYduig8c6b3opYfvuO57bdwDzmULBwAwKE1fJoyI2NPiffqeSBgLX1qE1gloAGAgaq3PR8TfNYx/qZRyTUu3eqvh860tXR9+hIAGAAbpfMd7/3ZL92h6Ar1+bt+hfnuwYVkENAAwSAcj4j8+8Nl7EfGNiHikpXs0PYGO8BSaARDQAMDAfOB47+MR8fsR8bFa692LWzyW7eE7rjsdzW/YsA+a1nmNHQAwaH8ZEdMR8eeL74gehLciYqbP555A0zpPoAGAgaq1nqm1PjHAeI7wJg6GSEADAKtB0z7oTXP7DukdWmULBwCwGvxvRMSpkye6z//TgY+dO3j+Hw+8tffOl9+JiP+stb43ktWxqghoAGA1OBkRvVdfOtT9zt88vqv25ud7tdervfn52uv9+uLP3BoRr45wjawSjvIGAFaFuX2HOnvv/NRlEXGg4UdurbUKaJZNQAMAq0YpZUcIaAbMpnoAYORKKTeUUv6qlPKJUa8FPoo90ADASJRSJmJhX/KeiPjM4seHI+JLI1sULIEn0ADA0JVSfjUivhcRfx0/jOeIiDtKKf0ORIEVQ0ADAKOwKSJ+ss/n68MTaFY4AQ0AjMK3ovmVcr9ZStkwzMVAhoAGAIZu8UCTP2oYz0TEF4e4HEgR0ADAqDweCweg9DNXStEprEj+YAIAI1FrPRERX28YfyIibslec3rb5ROdzsTERHdysju1fv3kuo0bpzZsnl63YfP0shYL5xDQAMAo/WFEzDfM9mQvdvW1P7d1asOmzZPrNmzsTk6tm+h2JzudTqd0Op3pbZdPLG+psEBAAwAjU2v9QUT8RcP4Z0spn85c7+gr33+naXb5T3xyY+Za0ERAAwCj9vB5ZnOZC73+g+dP14bZlm2Xe7MHrRDQAMBI1Vr/PSIONow/X0r5+FKvNX/mvYher9dvtnHLJQKaVghoAGAl+GrD552IeCBzoVr7B/S6DZsFNK0Q0ADASvD3sXC0dz93l1IuWeqFer3aN6An128U0LRCQAMAI1dr7UXzXuhNEXHv0q/V/wl0d2q9gKYVAhoAWCn2R8ThhtlvlVLWLeUijQHdndowt+9QudDFwf8T0ADAilBrPR0RjzaML4uILyzpOg1fIoyFkw2XFOFwPgIaAFhJ/iQiTjXM5kopH/kEuTGgF9jGwbIJaABgxai1HouIbzSMfzoifnFJV2l6GXTE+gtZF5xLQAMAK80jEdGUwB91vPfZiHitOzl1ZHLdhuPn/tWdWjcfnkDTgu6oFwAAcK5a68ullG9HxG19xp8tpXyq1vrdht89HBG/PLfv0A0Rsb3Pj3gCzbJ5Ag0ArET9DlY5Gwtv6ji5hN8/3fC5gGbZBDQAsOLUWr8TEf+y+I9vRcRXIuLjtdYv1lpfWsIlmgLaFg6WzRYOAGCl+r2ImI2Ir9da30r+btObPDyBZtkENADQirLwnuVtLV7yu4t/rSulzGR+8TO/8qVNV/7Up7d8aFDr5sl1N8ycfa/pAfVQHFs8eZExVc73nhcAgKVajNymkwT5oUtrrUdGvQgunD3QAACQIKABACBBQAMAQIIvEQIAA/PCCy/E9u39zjMZjPn5+ThypP/24pmZmZiYmBjaWiIijh49Gtdcc81Q78ngCWgAYGC2b98eMzOpF2ikPPvss/HMM8/Egw8+GBERr776atx55519f/bAgQNx2WWXDWwtrB22cAAAY6XX68XTTz8dN954Y1x//fXx0EMPxXPPPTfqZbGGCGgAYCz0er148sknY3Z2Nm655ZY4ePDg+7O9e/eObmGsOQIaABgLpZR4/PHH48UXX/zQbP/+/fHaa6+NYFWsRQIaABgLpZTYs2dP39mZM2fi0UcfHfKKWKsENAAwNm677ba46qqr+s4ee+yxePvtt4e7INYkAQ0AjI1utxsPPPBA39nx48fjqaeeGvKKWIsENAAwVu6666646KKL+s6eeOKJqLUOeUWsNQIaABgr09PTsXv37r6zV155JU6cODHkFbHWCGgAYOzcd999MTk52Xf2+uuvD3k1rDUCGgAYOzt27IidO3f2nb3zzjtx8uTJIa+ItURAAwBjaW5urnHmKTSDJKABgLE0Ozsbn/vc5/rO3nzzzTh9+vSQV8RaIaABgLHVdLBKRMQbb7wxxJWwlghoAGBs3XTTTXHttdf2nR07dizOnj075BWxFghoAGBsne9471prHD58eMgrYi0Q0ADAWLv99tvjiiuu6Ds7cuRI9Hq9Ia+I1U5AAwBjbXJysvF477Nnz8axY8eGvCJWOwENAIy9e+65J6anp/vOfJmQtgloAGDsbd26NXbt2tV39u677zrem1YJaABgVbj//vtjYmKi78zBKrRJQAMAq8KVV14Zt956a9/Z1NRUnDp1asgrYrUS0ADAqrF79+73/77T6cSll14as7OzcfXVV8eGDRtGuDJWk+6oFwAA0Jabb745du7cGbOzs3HvvffGxRdfPOolsQoJaABgVdm/f/+ol8AqZwsHAAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkNAd9QIAgNXr6NGjo17CSK31f//VqtRaR70GAGAVKKXMRMThUa9jDFxaaz0y6kVw4WzhAACABAENAAAJAhoAABLsgQYAWlFK6UTEtlGvYwwcq7X2Rr0ILpyABgCABFs4AAAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAASBDQAACQIaAAASBDQAACQIKABACBBQAMAQIKABgCABAENAAAJAhoAABIENAAAJAhoAABIENAAAJAgoAEAIEFAAwBAgoAGAIAEAQ0AAAkCGgAAEgQ0AAAkCGgAAEgQ0AAAkCCgAQAgQUADAECCgAYAgAQBDQAACQIaAAAS/g/swMdxavuadwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 750x750 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "points = [[0.20066186907158112, -0.032855948849496976], [0.3703214095102396, -0.1323611333831328], [0.48125568589863593, -0.29269832005189755], [0.5339926472234708, -0.48612209986502347], [0.5527041638633724, -0.6915413640867977]]\n",
    "points = [(0, 0)] + points\n",
    "points = np.stack(points, 0)\n",
    "points = points[:,::-1]\n",
    "target_index = 2\n",
    "\n",
    "c, r = ls_circle(points)\n",
    "arc = make_arc(c, r, start=(np.pi + 1.0 * np.pi), end=2 * np.pi + 0.62 * np.pi)\n",
    "\n",
    "# Perturb after to emphasize projection.\n",
    "points[target_index] += [0, 0.1]\n",
    "target = points[target_index]\n",
    "closest = project_point_to_circle(target, c, r)\n",
    "\n",
    "plt.figure(figsize=(2.5, 2.5), dpi=300)\n",
    "\n",
    "ax = plt.gca()\n",
    "\n",
    "x = Line2D([0, 0], [0, 0.5], color='k')\n",
    "y = Line2D([0, 0.5], [0, 0], color='k')\n",
    "\n",
    "project_to_arc = Line2D([target[0], closest[0]], [target[1], closest[1]], linestyle='-', color=colors['scarletred1'], alpha=0.5)\n",
    "to_target = Line2D([0, closest[0]], [0, closest[1]], linestyle='--', color='k', zorder=6.0)\n",
    "up = Line2D([0, 0], [0, np.linalg.norm(closest)], linestyle='--', color='k', zorder=6.0, alpha=0.8)\n",
    "angle_arc = get_angle_plot(x, to_target, offset=0.5, color='k')\n",
    "\n",
    "ax.plot(arc[:,0], arc[:,1], alpha=0.3, linestyle='--', color=colors['skyblue1'], zorder=0.0)\n",
    "\n",
    "\n",
    "ax.add_line(to_target)\n",
    "ax.add_line(project_to_arc)\n",
    "ax.add_line(up)\n",
    "ax.add_patch(angle_arc)\n",
    "\n",
    "\n",
    "# Car\n",
    "W = 0.1\n",
    "H = W * 1.9\n",
    "car = Rectangle(\n",
    "        (0 - W / 2, 0.0 - H * 0.8), W, H,\n",
    "        fill=True, edgecolor='k', facecolor='w', zorder=5.0)\n",
    "ax.add_patch(car)\n",
    "\n",
    "\n",
    "ax.text(-0.09, 0.32, r'$s^*$', fontsize=12)\n",
    "\n",
    "#for i in range(1, len(points)):\n",
    "#    ax.text(points[i,0] - 0.06, points[i,1], r'$p_%s$' % i, fontsize=8)\n",
    "ax.scatter(closest[0], closest[1], color=colors['skyblue1'], zorder=8, s=40, edgecolor='k')\n",
    "ax.scatter(points[1:,0], points[1:,1], color=colors['scarletred1'], zorder=8, s=40, edgecolor='k')\n",
    "\n",
    "#plt.plot(points[0,0], points[0,1], 'k.')\n",
    "\n",
    "BOUNDS = 0.7\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.axis('off')\n",
    "ax.set_facecolor((1.0, 1.0, 1.0))\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.gca().axis('equal')\n",
    "plt.xlim(-0.79, BOUNDS / 4)\n",
    "plt.ylim(-0.02, 0.25)\n",
    "\n",
    "plt.savefig('controller.pdf', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.283185307179586 8.230972752405258\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAALQCAYAAAC5V0ecAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAuIwAALiMBeKU/dgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+M3Odh5/fPl6SXXFrWDiMpMW0ppqREP76ORessN+xZl7vAvgLx4VCgdwYqJSlS1LL+iAM0LpwDAtgFcvkjTVynh1PRKjL6Tw/SH8b1egEatD2rvWujmjlLkWnHX1mOFVGmJCq2HM5a0i5JcfXtHzO0doc75DzL+bW7rxcwEHe+s995ZErWm88+8zxV27YBAABGs2fWAwAAgO1EQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBg36wHsNtUVbUvyT39x4eS3JhkIcn5JC8meSrJk0mebNv2wqzGCQDA5qq2bWc9hl2hqqp3J3kgyYNJ3jvCt7yU5OEkj7Rt+8okxwYAwOgE9IRVVbU3yWeS/E6SA1u4xdkkn0/yxbZt18Y5NgAAygnoCaqq6uYkjyX5+THc7niS+9u2fX4M9wIAYIsE9IRUVVUn+UqSw5tdv31hIXcvLuaO/ftzsKqy0rb59rlzeXp1Nc+ePz/stqeTfKxt22ZCwwYA4AoE9AT0Z56fyEA8V0nuX1rKA4cO5eji4tDvP7G6mkfOnMmjy8vZ5HfndJKPmIkGAJgNAT1m/TXPT2Rg2catCwt56PDhHDt4cOR7HV9ZyadPn85zl85IH09yrzXRAADTZx/o8ftMBuL5w4uLefzIkaJ4TpJjBw/m8SNH8uFLZ6uPJfnNqxolAABbYgZ6jPpb1T2fdbtt3LqwkMePHMnS3r1bvu/y2lo+evLk4Ez02SQ32+IOAGC6zECP1wNZF89VkocOH76qeE6Spb1789Dhw6k2Pn0gySev6sYAABQT0GPSP2HwwfXP3b+0VLxsY5hjBw/mvqWlwacf7L8vAABTIqDH554MnDD4wKFDY32DT116vxvTOw4cAIApMXs5Pves/+L2hYXLblW3FUcXF3PbwkK+s24t9HV79/6dbl0/1WmaC2N9M+ZG/6cM9/QfH0rvD04LSc4neTHJU0meTPJk27b+OQCACRPQ47NhJvjuMcfz+vuuD+ifXVj45STXdOv6XJLX+483Lvdrsb099D+U+kB6S4Pee5mX/lr/ry9VVfVwkkd8uBQAJkdAj8+N67+4Y//+ibzJnQP3/eHa2rX9X+7vP6670j26df3vO03zJyO8ruo0jW1apqy/l/hnkvxO1n0odQTv7X/Pb1dV9fkkX7RXOACMn4Aen4X1XxysqmGvuyqLA/e90LZb2eLj7IivO9at67+bK8xoX/y1me2r1z/F8rEM7CVe6ECS30/yn1RVdb9TKwFgvAT0+GzYpHllQvtrrw7cd19VbWWG8fURX3dNejF2IKPNbJ/NCEtIIrY3VVVVneQrGTgC/qLbFxZy9+Ji7ti/PwerKittm2+fO5enV1fz7KWnVSa9A3eeqKrqY23bNhMcOgDsKgJ6fF5c/8W3z52byJs8M3Df6/bu/dEWbvPGiK97Z+F9L8b29Vd6Ybeuj3ea5n8f4XW7YhlJf+b5kniu0tsO8YFDhy77odQTq6t55MyZPLq8nIH/sQ4n+UpVVR8xEw0A4yGgx+epvP1hrjy9ujqRNxm8b71//8tbuE3JDPSkjPonjP+wW9e/kNE/ILnt1vz21zw/loF4vnVhIQ8dPjzSXuJHFxfz0OJifqXTyadPnx48tfJwkkerqrrXmmgAuHoCenyeXP/Fs+fP58Tq6li3svv66uqGHTiS5Jvnzv2zJN9Jb7b4mv5js1+vXys96gz0JAN6K8tIRpnZXs3oa7bnJSY/k4E1zx9eXMyXb7qp+BTLYwcP5vEjR/KJU6fytY1/2DqW5DeTfOFqBwsAu52AHp8nk7yUdduNPXLmTB4aY0A/cubM4FMvfm119f+60nribl1X6QXoxZjujviWpUs4SkxqFnyx/xgltv+/TtP8nyO8bmLLSPpb1f3O+uduXVjYUjxftLR3b75800356MmTgzPR/7Sqqn9hizsAuDoCekzatr3Q34P3xzH06PJyfqXTGctx3l9dWcljy8uDTz88ysEZ/fhb7T9eLXjbf5m3Z7KHzWxv9TTLeZgFH3UZyd/u1vVHMsKHI1M+s/1A1m1VVyV56PDhLcfzRUt79+ahw4fz8RdeWL8m+kCSTyb53au6OQDscgJ6vB5J8tvpB1Gb5NOnT+fxI0euKoiW19byG6dPD3447GySL235piPoNM3Jy13vz2wv5u2YvtIykvWxPeoM9CRnwUeN+HclOdh/3HClF3freiUjLCP5rVdeOZfeISk/dv/S0lj+wJX0lnPct7SURzf+wevBqqp+z4mFALB1AnqM2rZ9pX+Axe9ffO658+fziVOntvwj+eW1tXzi1KnBH8Unyedm/aP4/sz2Sv/xg8u9dpPYHnX3kHlYh10a8SPF9j2Lizc+cubMhhMGHzh0qPCtLu9Thw4NBvSN6Z2a+WdjfSMA2EUE9Ph9Mck/yroPhX1tdTUfPXly5B0VLjq+srLZjgpJcjzJH45jsNNSEtsDvpzeDPD6mezB2e2tLiOZ6W4kx1dW3rP+69sXFsb6odOktzvHbQsLgx8+vScCGgC2TECPWdu2a1VV3Zfkiazbluy58+fz8RdeyH1LS/nUCHv6/tGZM3ns0j19k+R0kvt3y3ZkBctIRllCMhjbM12H3Zw7t2HburvHHM/r7zsQ0B+ayBsBwC4hoCegbdvnq6r6WAYOxmjT+2Dho8vLua1/qtyd+/dnsaqy2rZ5pn+q3OBWdeucTvIxB2K8bWBm+7IGYrtkGclE1mH/zdrateu/vmP//km8Te689L43TuSNAGCXENAT0rZtU1XVR5I8mt4evBt85/z5y4XyZo6nN/MsnrdoILa/X/CtX87ldyPZ0jKSC227YVH8waoqvcVIFi+972RKHQB2CQE9Qf2Z6HvTOyjjd7Juu7ICZ5N8Lskf7pZlG/Om0zSX/UNLf2b7YK68hOTir6sk2VdVG34/V9rJnFi+OnDfn9y7953dur4lyfeutIc4AHApAT1h/ej9g6qq/uf09uB9MKP9CP3FJA8n+dKsd9vg8voz22/0H5ed2e7W9Z70l5G8urb2/iS3XLz27XOjbktd5pmB+96ysLA3yX+W5M1uXb+Q5E+vtNYcAHibgJ6SfgT/blVVv5feh7ju6f/1xvR+pH4uvWh+Kr1TDZ+yV+/O02mat9KP7TNV9e+S/OOL157eePT22Azet96//+X+L9+R5GeSfHUibwwAO5SAnrJ+FP9ZbCNG7w9KP/bs+fM5sbo61q3svr7Jh1KPHTz48rovLyT53tjeEAB2AQENs/NkkpeS/PgwlUfOnMlDYwzoR86c2fD1tXv2/Ogfvutdp9c99b1O07x5pfv013n/4/R+SvJckh/0l64AwK4joGFG2ra9UFXVw+l9wDRJb5vDX+l0xnKc91dXVvLYxlMI8/evueapA3v2vLXuqedGvN1PJXl//5Ekr3Xr+rkkf5XkrzpNM+qhNACw7QlomK1Hkvx2+ju0tEk+ffp0Hj9yZEtHv1+0vLaW3zh9esNBPFVy7r84dOixJNfl7R1hRg3oWwa+fleSD/Yf6db1K/17PRe7ewCww1XthLbOAkZTVdVnk/z++uc+vLiYL99005YienltLZ84dSpfu/RDiZ9t2/YL/Z1A3pPkSJInRlmK0a3rX01y64hDeDPJC3k7qC33AGBHEdAwY1VV7U3v6PefX//8rQsLeejw4aLlHMdXVvLp06fz3KWH9BxPcu9W9hLv1vU7kvyTbP0nVq+lF9KWewCwIwhomANVVd2cXkQf3vB8kvuWlvKpQ4cuuzvHidXV/NGZM3lseTmb/Bt9OslHtnqKZbeub03yq1v53iEs9wBgWxPQMCeqqqqTfCUDEX3RbQsLuXtxMXfu35/Fqspq2+aZc+fy9CZb1a1zOsnH2rZttjqubl0fSnJXeks4bswWji2/jDeT/C+dpnlmjPcEgIkS0DBH+jPRjyY5NobbHU9y/1ZnnjfTrev96a2dvrX/uG4Mt32o0zSvjuE+ADAVAhrmTH9N9GfS297uwBVevpmzST6X5A+3sua5RH92+pb0YvqWlI93Ocl/N+IHGfd2mmaifz8AMAoBDXOqqqp3J/lkkgfTWzpxJS8meTjJl/pHx0/Vut09Ls5Oj7Lc4887TfPHI9x7f3p/qDiV3ocRn0vyfbt7ADALAhrmXFVV+5J8KMk9/b/emGR/knPpRfNT6Z1q+FT/qPi5MOJyjy93muZbI9zr9iT3DTz9Wt6Oabt7ADA1AhqYim5dd/J2TN+S3h8C/qDTNCsjfO8vZWCbv03Y3QOAqRDQwNT1l3vc0Gmavx7x9Z9Ocn3BW1zIxsNcLPcAYGwENDDX+jPX/+VV3ub1vB3TlnsAcFUENDDXunV9d5L/eMy3vbjc46/SW+7x5pjvD8AOttWjeQGm5USSH6Rsd48reXf/8ZEkjyV59irvB8AuYgYa2FbGfJjLW0n+m07TnBvD0ADYJQQ0sK1tsrtHyWEu3+s0zf80wntUSd6X5MVZ7+7R39bwnmzc1nAhyfls3NbwyXna1hBgJxHQwI6x7jCXi6cj3pTLL/f4vztN8+9GuO91SX4jM9zdo3+wzgPpHazz3hG+5aX0DtZ5ZBYH6wDsZAIa2LFGWO7xpU7TvDjCff6DJB/f5NLEd/cY09Hun0/yxUkf7Q6wWwhoYNcYWO5xOMk/7zTNWyN8331Jbh/hLcZ6mEtVVTen9yHHKx0iM4rjSe5v2/b5MdwLYFcT0MCu1K3rapTlF9263pvkn6S3zrjEVS33qKqqTvKV9EL/ErcvLOTuxcXcsX9/DlZVVto23z53Lk+vrubZ8+eH3fZ0ko+1bdsU/Z0AsIGABriMbl2/L8l/PoZbjbzcoz/z/EQG4rlKcv/SUh44dChHFxeHvtGJ1dU8cuZMHl1ezib/D386yUfMRANsnYAGuIxuXf9ckn+QZHixbs0r6R3k8rVO05y5+GR/zfMTGVi2cevCQh46fDjHDh4c+Q2Or6zk06dP57lLZ6SPJ7nXmmiArRHQAFfQ393jcN5eP32l3T1K/I+dpvnxLhlVVX02ye+vf8GHFxfz5ZtuytLevcU3X15byydOncrXVlcHL322bdsvbGXAlLP9IOwsAhqg0BgPc3kjyRcuro3ub1X3fNbttnHrwkIeP3JkS/F80fLaWj568uTgTPTZJDfb4m6ybD8IO5OABrhK63b3uKX/GHW5xzc7TfMvL35RVdXn0tuurvd1kj953/uKlm0Mc3xlJR9/4YXBNdGfa9v2d6/65lzC9oOwswlogDEqXO7xv3aa5uvJj3/EfzLrZil/eWkpD73nPWMb26+//HIeXV5e/9SL6c1CWzIwRrYfhJ1PQANM0BWWe/y3naZ5LUmqqjqW5Kvrv/ffHjly2d02Sp1YXc3fO3ly8Oljbdv+2djeZJez/SDsDvtmPQCAnazTNOeSPNt/rF/uccPFeO67Z/333b6wMNZ4TpKji4u5bWEh31kXajfs3fsLSQT0GPRnni+J5zFsP3g4yVeqqrL9IMwJAQ0wRZ2m6aa348KgD63/4u4xx/P6+64P6NW2vXuU7+vW9a8muTbJav9xdsiv1399ttM0u2L9bn/N82MZiOeS7QePLi7mocXF/Eqns9n2g4eTPFpVle0HYQ4IaID5cOP6L+7Yv38ib3LnwH3PvvXWT434rdcnWSp9v25dn8vwwL7462cHZuO3o89kYM3zVrcfPHbwYB4/cmSz7QePJfnNJLYfhBkT0ADzYcNR4QeraiJvsjhw37XkHSN+61Z2kkiS/f3H5eL7+0muGNDduv6FJO/KFWa+O03z5hbHuiX9rep+Z/1zty4sbHnv7iRZ2rs3X77pps22H/ynVVX9C1vcwWwJaID5sKGSVib0Ae/Vgfu2vei8rG5d700vgiflimPoq5O8+0ov6tb1hVx+1nuzWfDlqwjvB7LuDxhVkocOH76qvbuTXkQ/dPjw4PaDB5J8MontB2GGBDTAfHhx/RffPnduIm/yzKX3fXGz1w3Y6uzzqC45JnGIUReG70tyTf8xqsfS/6Dn5XTr+tb0/vdYTbJ64uzZN9M7JOXH7l9aGsve3UlvOcd9S0uD2w8+WFXV79l+EGZHQAPMh6eS/NrFL56+9Ojtsdjkvpt9oHHQdgvoSY7h76S3LWGS5Nlz527MwAmDDxw6NL5RJfnUoUODAX1jeh86tXsKzIiABpgPT67/4tnz53NidXWsW9l9fXV1ww4cm73vECtJ/lV6AbuYXlAvDvl62KExw7zZaZorzqT2l5EsXOl1V2HUgN7wh4njKysbTrqZ1vaD6W17KKBhRgQ0wHx4MslLWTeb+ciZM3lojDH2yJkzg0+9mBFmoDtNs5rkxJVe163rKr0PJQ6L683Ce9Qt2SY9Cz7qOuwNvyHNuXMbtq2b1vaDGdj2EJguAQ0wB9q2vVBV1cNZt5vDo8vL+ZVOZyzrab+6spLHNi4DSJKHx7mOttM0bXofhjyf5JI3G4Ons3mcj7qTyOVsaRnJ36ytXbv+62ltP5iBbQ+B6RLQAPPjkSS/nf5sa5vk06dP5/EjR65qR4fltbX8xunTg6fbnU3ypS3fdMo6TfNGkn+92bVuXe/LlWe6h13bk6tYRnKhbTf8xkxr+8FMdlcU4AoENMCcaNv2laqqPp/k9y8+99z58/nEqVNb3lN4eW0tnzh1anAv4ST53E7ZS7gfv6/3HyNbt+Rk1OUhC0nO5O0Qz76q2rAEZVrbDyaZzDYtwEgENMB8+WKSf5R1p9p9bXU1Hz15cuQjoS86vrKy2ZHQSXI8yR+OY7Db2cCSk1Fev5rknyVJt673JDnw4ptv/mySWy6+Zs62HwQmpPTT0gBMUNu2a0nuS3J6/fPPnT+fj7/wQn795Zdz4gpb3J1YXc2vv/xyPv7CC5vF8+kk9/ffhy3qNM1bnaZZWW3br65/fs62HwQmpGon9OMmALauqqo6yVeSHN7s+m0LC7l7cTF37t+fxarKatvmmXPn8vTmW9VddDrJx9q2bSY07F2nqqpjSTZE9L89cmTs2w/+4smTg08fa9vWNnYwIwIaYE5VVXVzkkeTHBvD7Y6nN/P8/BjuRV9VVfuSnMy67Qd/eWkpD73nPUO/p9Svv/zy4EEqLya52UmEMDuWcADMqX7s3pvktzL6PsWDzib5bJJ7xfP49SP24fXPPbq8nOMrK2O5/zS2HwTKmYEG2Aaqqnp3kk8meTCj7QH8Ynph96WdstvGvOr/3jyfdbt53LqwMJbtBz968uTgOvaz6c0++z2FGRLQANtIf8nAh9I7yvlD6cX0/vS2Nbt4suCTSZ4ySzk9VVV9Nuu2H0ySDy8uXvX2g1+79MODn23b9gtbHykwDgIaAK5SVVV7kzyRddsPJr2Z6DFvP3ivHVRg9gQ0AIxB/0OfT2Rg55QqyX1LS/nUoUOX3Z3jxOpq/ujMmTy2vDx4amTS20HlI9axw3wQ0AAwJrYfhN1BQAPAGNl+EHY+29gBwBiNY/vBfcmFpT17PhfbD8JcMgMNABNSuv3gtXv2/OjvX3PNU//Vddf9+Z0HDvzrTtP8vxMfJFBMQAPAhF1u+8EPHjhwzd86cODMsYMHX/6H73rX6QN79rzV/7ZXk/z3nabxH2qYM/tmPQAA2On6e3L/Wf+xQbeu/3aS/2iTb7s+yXuSvDTZ0QGlrIEGgNn6ZrLZznVJkqPTHAgwGgENADPUaZrXkjw35PIHunXtp8UwZwQ0AMzeiSHPLyb52WkOBLgyAQ0As/ftJOeGXLOMA+aMgAaAGes0zZtJvjXk8s926/rgNMcDXJ6ABoD5MGwZx94kPzfNgQCXJ6ABYD58L8mZIdc+OM2BAJcnoAFgDvQPTPnGkMvv6db1DdMcDzCcgAaA+TFsGUfiw4QwNwQ0AMyJTtP8TXpLOTZzV7eu/Xcb5oB/EQFgvgybhb42yc3THAiwOQENAPPlW0kuDLlmGQfMAQENAHOk0zRn0ztYZTN3dut6/zTHA1xKQAPA/Bm2jOMdSe6c5kCASwloAJg/zyV5fcg1yzhgxgQ0AMyZTtO8leSbQy7f3K3rzjTHA2wkoAFgPl1uT+i7pjYK4BICGgDmUKdpXknyysDTbyT5aoZ/yBCYgn2zHgAAMNSJJNcnebb/6+/2l3cAMySgAWB+/XmSr3eaZnXWAwHeVrVtO+sxAADAtmENNAAAFBDQAABQQEADAEABAQ0AAAUENAAAFLCNHQBsU9263pPkSJIznaY5M+PhwK4hoAFgm+nW9Q1JjqZ3pPe1SZ5I8m9mOijYRewDDQDbRLeu70jyd5K8d+DSa0n+0CmFMB3WQAPA9nEwl8ZzkrwryS1THgvsWgIaALaPbyW5MOTa0WkOBHYzAQ0A20Snac4leWbI5Tu7db1/muOB3UpAA8D2cmLI8/uSvH+aA4HdSkADwPbyV+l9aHAzlnHAFAhoANhG+jttfGPI5fd16/rQNMcDu5GABoDtZ9gyjqS3NzQwQQIaALaZTtN8P8npIZePduu6muZ4YLcR0ACwPQ2bhf6JJDdNcyCw2whoANievplk2MmDPkwIEySgAWAb6jTNG0n+csjl93fr+h3THA/sJgIaALavYcs4DiS5fZoDgd1EQAPA9vWdJKtDrlnGARMioAFgm+o0zYUkfzHk8q3dur5mmuOB3UJAA8D2NmwZx54kH5jmQGC3ENAAsL29lOSHQ659cJoDgd1CQAPANtZpmjbDZ6F/qlvX757meGA3ENAAsP1d7mhvHyaEMRPQALDNdZpmOcnzQy7f1a1r/72HMfIvFADsDJvNQp9O8v8k2TvlscCOtm/WAwAAxuKZJP8gydkk30hyotM035/tkGBnqtq2nfUYAIAx6Nb1TyX5Qadp3pr1WGAnE9AAAFDAGmgAACggoAEAoICABgCAAgIaAAAKCGgAACggoAFgl+jW9f5uXb9v1uOA7c42dgCwg/WP8b4lydEkdyZZS/KFTtO8OdOBwTbmJEIA2IG6dX0gyS8k+UCSd627tC/JHUm+OYtxwU5gCQcA7ExvJrk7G+P5oqNTHgvsKAIaAHagTtOsZfgs863dut4srIERCGgA2LlODHm+SnLXNAcCO4mABoCd6+Ukrw65drRb19U0BwM7hYAGgB2q0zRtkq8PufyTSd49xeHAjiGgAWBn+0aSYXvW+jAhbIGABoAdrNM0P0ry/JDLH+jW9d5pjgd2AgENADvfsA8TvjPJz0xzILATCGgA2PmeSXJ+yDXLOKCQgAaAHa7TNOeTNEMu396t68Vpjge2OwENALvDsGUce5P83DQHAtudgAaA3eFkkuUh1yzjgAICGgB2gf6e0N8YcvnGbl1fP83xwHYmoAFg9xi2jCNxtDeMTEADwC7RaZpXk7w45LKjvWFEAhoAdpdhs9BLSY5McRywbQloANhd/iLJ2pBrPkwIIxDQALCLdJpmNcmzQy7X3bpemOZ4YDsS0ACw+wxbxrGQ5M5pDgS2IwENALvPd5OsDDy3luTbSc5MfziwvVRt2856DADAlHXr+peS/HySl9Kbkf6LTtMMRjWwiX2zHgAAMBNfTfJkp2l+MOuBwHZjBhoAAApYAw0AAAUENAAAFBDQAABQQEADAEABAQ0AAAUENACwqW5d/0S3rg/Mehwwb+wDDQD8WD+Y35/kaJKfTvInSf79TAcFc0ZAAwDp1vWRJB9Ocns29sHRCGjYwBIOACBJbk5v5nlwcu293bq+YQbjgbkloAGAJDlxmWtHpzYK2AYENACQTtP8TZJTQy7f1a1rzQB9/mUAAC76+pDnr01yZIrjgLkmoAGAi76V5MKQa5ZxQJ+ABgCSJJ2mOZvk2SGX625d75/meGBeCWgAYL1hHyZ8R5I7pzkQmFcCGgBY77tJ3hhyzTIOiIAGANbpNM1bSb4x5PLN3bruTHM8MI8ENAAw6HJ7Qn9gaqOAOSWgAYANOk3zSpK/HnL5g926rqY5Hpg3AhoA2MywWejrkrx3mgOBeSOgAYDNfCNJO+SaDxOyqwloAOASnaZ5Pb0dOTbzc9263jfN8cA8EdAAwDDDlnEsJrltmgOBeSKgAYBhnk1ydsg1yzjYtQQ0ALCpTtO8meRbQy7/bLeu3znN8cC8ENAAwOUMW8axJ8nPTXMgMC8ENABwOaeSnBlyzTIOdiUBDQAM1WmaNpvPQr+R5Hvdut475SHBzNmCBgC4khNJ/l6StfQ+WHgiyXc7TbM2y0HBrFRtO2yPdACAnm5d10me7zTN6qzHArMmoAEAoIA10AAAUEBAAwBAAQENAAAFBDQAABQQ0AAAUEBAAwBXrVvXVbeul2Y9DpgG29gBAFvWrevr0zvS+670Dlr55/3TC2HHchIhAFCkW9f70wvmo0luHLh8Y5JTUx8UTJElHABAqX1JfimXxnOSfHDKY4GpE9AAQJFO07yR5LtDLr+/W9d+ws2OJqABgK34+pDnDyS5fZoDgWkT0ADAVnwnydkh145OcyAwbQIaACjWaZoLSf5iyOWf6db1NdMcD0yTgAYAturEkOf3JPnANAcC0ySgAYCtejHJD4dcs4yDHUtAAwBb0j8wZdgs9Lu7df1T0xwPTIuABgCuxjcuc80sNDuSgAYAtqzTNN0kJ4dcvqtb11qDHcc/1ADA1Rq2jOOaJLdOcyAwDQIaALhaTZI3h1yzjIMdR0ADAFel0zTnkjwz5PI4xOw2AAAREElEQVQd3bo+MM3xwKQJaABgHIYt49iX5P3THAhMmoAGAMbh+SQ/GnLNMg52FAENAFy1TtO8leFb2v10t65/YprjgUkS0ADAuAxbxpEkd01tFDBhAhoAGItO0/wgyctDLh/t1nU1zfHApAhoAGCchs1CH0ry09McCEyKgAYAxumbSd7a5PlX0tuRA7a9qm3bWY8BANhBunX9nya5I8nr6X2w8ESnaf56tqOC8fEnQQBg3P40yVNJnuvvzgE7ihloAAAoYA00AAAUENAAAFBAQAMAQAEBDQAABQQ0AAAUENAAwNR163p/t651CNuSfaABgKnoB/PNSY4muTPJl5N8Z6aDgi2wDzQAMFHdur4hvWi+K8m1Z996a88fv/bae/74Rz+q/rfXX9+T5MYkC0nOJ3kxvUNYnkzyZNu2F2Y1bhjGDDQAMGn3Jjn6zNmz13zhhz/8W//m9dfvee2tt951mdf/Wv+vL1VV9XCSR9q2fWXio4QRmYEGACbqL2+77db/+vvf/x++vLz8ixe2Nnl3Nsnnk3yxbdu1MQ8PigloAGBiqqq6OcljSX5+DLc7nuT+tm2fH8O9YMsENAAwEVVV1Um+kuTwZtdvX1jI3YuLuWP//hysqqy0bb597lyeXl3Ns+fPD7vt6SQfa9u2mdCw4YoENAAwdv2Z5ycyEM9VkvuXlvLAoUM5urg49PtPrK7mkTNn8ujycjYpldNJPmImmlkR0ADAWFVVtTe9eN6wbOPWhYU8dPhwjh08OPK9jq+s5NOnT+e5S2ekjye515poZsEG5gDAuH0mA/H84cXFPH7kSFE8J8mxgwfz+JEj+fCls9XHkvzmVY0StsgMNAAwNlVVvTvJ80kOXHzu1oWFPH7kSJb27t3yfZfX1vLRkycHZ6LPJrnZFndMmxloAGCcHsi6eK6SPHT48FXFc5Is7d2bhw4fTrXx6QNJPnlVN4YtENAAwFhUVbUvyYPrn7t/aal42cYwxw4ezH1LS4NPP9h/X5gaAQ0AjMs9Sd67/okHDh0a6xt86tL73ZjkQ2N9E7gCAQ0AjMs967+4fWHhslvVbcXRxcXctrBw2feFSRPQAMC4bJgJvnvM8XyZ+5qBZqoENAAwLjeu/+KO/fsn8iZ3XnrfGzd7HUyKgAYAxmXD2oqDVTXsdVdl8dL7TqbUYQgBDQCMy4ZNmlcmdNbE6qX3PTeRN4IhBDQAMC4vrv/i2+cm07XPXHrfFzd7HUyKgAYAxuWp9V88vbo6kTfZ5L5PbfY6mBQBDQCMy5Prv3j2/PmcGHNEf311Nd/ZeJz3Je8LkyagAYBxeTLJS+ufeOTMmbG+wSb3ezFmoJkyAQ0AjEXbtheSPLz+uUeXl3N8ZWUs9//qykoeW14efPrh/vvC1FTthD4hCwDsPlVVvTvJ80kOXHzu1oWFPH7kSJb27t3yfZfX1vLRkyfz3MblG2eT3Ny27StbvjFsgRloAGBs+jH7+fXPPXf+fD5x6lSW19a2dM/ltbV84tSpwXhOks+JZ2bBDDQAMFZVVe1N8kSSn1///K0LC3no8OEcO3hw5HsdX1nJp0+f3iyejye5t23brVU5XAUBDQCMXVVVN6cX0Yc3PJ/kvqWlfOrQoRxdXBz6/SdWV/NHZ87kseXlbFIqp5N8pG3b58c6aBiRgAYAJqKqqjrJVzIQ0RfdtrCQuxcXc+f+/Vmsqqy2bZ45dy5Pb75V3UWnk3ysbdtmQsOGKxLQAMDE9GeiH01ybAy3O57kfjPPzJoPEQIAE9OP3XuT/FZ6u2Zsxdkkn01vzbN4ZubMQAMAU9Hf4u6TSR5McuMI3/JievtKf8luG8wTAQ0ATFVVVfuSfCjJPf2/3phkf5JzeftkwSeTPOWQFOaRgAYA5kK3rg8muT7JDf3H652m+dPZjgoutW/WAwAAdrduXX8oyS8muWbg0ukkApq540OEAMCsXcil8Zwk13frupr2YOBKBDQAMGs/GPL8O5J0pjkQGIWABgBm7dXLXLthaqOAEQloAGCmOk1zPkl3yGUBzdwR0ADAPBi2jENAM3cENAAwD4Yt4xDQzB0BDQDMg6Ez0HbiYN4IaABgHgwL6IUk105zIHAlAhoAmAfDAjqxjIM5I6ABgJnrNM3ZJK8NuSygmSsCGgCYF8Nmoa+f6ijgCgQ0ADAvbGXHtiCgAYB5YScOtgUBDQDMi2EBvZjkndMcCFyOgAYA5oWdONgWBDQAMBc6TbOS5I0hlwU0c0NAAwDzxAcJmXsCGgCYJwKauSegAYB5IqCZewIaAJgnwwL6nd26PjjVkcAQAhoAmCd24mDu7Zv1AAAA1nkjyWqSc+nF9PrHX89wXPBjVdu2sx4DAMCPdet6X6dpLsx6HDCMgAYAgALWQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABQQ0ALBtdOtauzBz9oEGAOZOP5SvS+/47vWP65L8Qadpzs1weOxyjvIGAObRUpJfH3Lt+iQvTXEssIEfgwAA86ib5M0h126Y5kBgkIAGAOZOp2naJK8OuXz9NMcCgwQ0ADCvfjDkeTPQzJSABgDmlYBmLgloAGBeDQvoQ926fsdURwLrCGgAYF4NC+gqve3sYCYENAAwr84kuTDkmmUczIyABgDmUqdp3krywyGXBTQzI6ABgHnmg4TMHQENAMwzAc3cEdAAwDwbFtA/0a3rfVMdCfQJaABgng0L6D1JfmKaA4GLBDQAMM/+JslbQ65ZxsFMCGgAYG51mmYtduJgzghoAGDe+SAhc0VAAwDzTkAzVwQ0ADDvhgX0dd261jJMnX/oAIB5Nyyg98ZOHMyAgAYA5t0Pk7RDrlnGwdQJaABgrnWa5kJ629ltRkAzdQIaANgONlvG8WaSd0x7IOAITABgO/huktX0QvriY7nTNMOWdsDEVG3rnzsAABiVJRwAAFBAQAMAQAEBDQAABQQ0AAAUENAAAFBAQAMAQAEBDQAABRykAgBsW926XkzvOO/lTtMsz3o87A4CGgDYNrp1fVeSG9OL5huSXNO/9H8k+eqsxsXuIqABgO3kaJJbN3n+hmkPhN3LGmgAYDv5wZDnBTRTI6ABgO1kaEB367qa6kjYtQQ0ALCdDAvoA3l7PTRMlIAGALaTVy9zzTIOpkJAAwDbRqdpVpK8PuSygGYqBDQAsN34ICEzJaABgO1GQDNTAhoA2G4ENDMloAGA7WZYQB/s1vU7pzoSdiUBDQBsN8MCOjELzRQIaABgu1lJsjrkmoBm4gQ0ALCtdJqmjXXQzJCABgC2IwHNzAhoAGA7EtDMjIAGALajYQF9TbeuF6c6EnYdAQ0AbEd24mBmBDQAsB29luTckGsCmokS0ADAtmMnDmZJQAMA29WwgL5+qqNg1xHQAMB2ZQaamRDQAMB2NRjQ3SR/meRb3bquZjAedol9sx4AAMAWvZTkX6UX0q92mub8jMfDLlG1bTvrMQAAwLZhCQcAABQQ0AAAUMAaaABgLKqq2pPkulmPYxv4Ydu2b816EGydgAYAxuW6JN+f9SC2gZ/M5Y8iZ85ZwgEAAAUENAAAFBDQAABQwBpoAGBimqbJ9ddfP+thzMyrr76auq5nPQzGTEADABNz/fXX54Ybbpj1MGCsLOEAAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKLBv1gMAAHauV199ddZDmKnd/ve/U1Vt2856DADADlBV1Q1Jvj/rcWwDP9m27Q9mPQi2zhIOAAAoIKABAKCAgAYAgALWQAMAY1FV1Z4k1816HNvAD9u2fWvWg2DrBDQAABSwhAMAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAoIaAAAKCCgAQCggIAGAIACAhoAAAr8/66NLwTMp+M2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 750x750 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "points = [[0.20066186907158112, -0.032855948849496976], [0.3703214095102396, -0.1323611333831328], [0.48125568589863593, -0.29269832005189755], [0.5339926472234708, -0.48612209986502347], [0.5527041638633724, -0.6915413640867977]]\n",
    "points = [(0, 0)] + points\n",
    "points = np.stack(points, 0)\n",
    "points = points[:,::-1]\n",
    "target_index = 2\n",
    "\n",
    "c, r = ls_circle(points)\n",
    "arc = make_arc(c, r, start=(np.pi + 1.0 * np.pi), end=2 * np.pi + 0.62 * np.pi)\n",
    "\n",
    "# Perturb after to emphasize projection.\n",
    "points[target_index] += [0, 0.1]\n",
    "target = points[target_index]\n",
    "closest = project_point_to_circle(target, c, r)\n",
    "\n",
    "plt.clf()\n",
    "plt.figure(figsize=(2.5, 2.5), dpi=300)\n",
    "\n",
    "ax = plt.gca()\n",
    "\n",
    "x = Line2D([0, 0], [0, 0.5], color='k')\n",
    "y = Line2D([0, 0.5], [0, 0], color='k')\n",
    "\n",
    "project_to_arc = Line2D([target[0], closest[0]], [target[1], closest[1]], linestyle='-', color=colors['scarletred1'], alpha=0.5)\n",
    "to_target = Line2D([0, closest[0]], [0, closest[1]], linestyle='--', color='k', zorder=6.0)\n",
    "up = Line2D([0, 0], [0, np.linalg.norm(closest)], linestyle='--', color='k', zorder=6.0, alpha=0.8)\n",
    "angle_arc = get_angle_plot(x, to_target, offset=0.5, color='k')\n",
    "\n",
    "# ax.plot(arc[:,0], arc[:,1], alpha=0.3, linestyle='--', color=colors['skyblue1'], zorder=0.0)\n",
    "\n",
    "\n",
    "# ax.add_line(to_target)\n",
    "# ax.add_line(project_to_arc)\n",
    "# ax.add_line(up)\n",
    "# ax.add_patch(angle_arc)\n",
    "\n",
    "\n",
    "# Car\n",
    "W = 0.1\n",
    "H = W * 1.9\n",
    "car = Rectangle(\n",
    "        (0 - W / 2, 0.0 - H * 0.8), W, H,\n",
    "        fill=True, edgecolor='k', facecolor='w', zorder=5.0)\n",
    "ax.add_patch(car)\n",
    "\n",
    "\n",
    "# ax.text(-0.09, 0.32, r'$s^*$', fontsize=12)\n",
    "\n",
    "#for i in range(1, len(points)):\n",
    "#    ax.text(points[i,0] - 0.06, points[i,1], r'$p_%s$' % i, fontsize=8)\n",
    "# ax.scatter(closest[0], closest[1], color=colors['skyblue1'], zorder=8, s=40, edgecolor='k')\n",
    "ax.scatter(points[1:,0], points[1:,1], color=colors['scarletred1'], zorder=8, s=40, edgecolor='k')\n",
    "ax.plot(points[:,0], points[:,1], '--', color=colors['scarletred1'], alpha=0.5)\n",
    "\n",
    "for i in range(len(points) - 1):\n",
    "    x = (points[i,0] + points[i+1,0])/2\n",
    "    y = (points[i,1] + points[i+1,1])/2\n",
    "    #ax.text(x, y, r'$s_%s$' % (i+1))\n",
    "    \n",
    "#plt.plot(points[0,0], points[0,1], 'k.')\n",
    "\n",
    "BOUNDS = 0.7\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.axis('off')\n",
    "ax.set_facecolor((1.0, 1.0, 1.0))\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.gca().axis('equal')\n",
    "plt.xlim(-0.79, BOUNDS / 4)\n",
    "plt.ylim(-0.02, 0.25)\n",
    "\n",
    "plt.savefig('speed.pdf', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
